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Abstract 


The  High  Reso]i-tion  Radar  Branch  of  the  Rome  Air  Development  Center  has 
developed  a tactical  target  identification  (TTI)  pulsed -Doppler  radar 
system  which  generates  two-dimensional  "Images"  of  aircraft.  The  signal 
processing  technique  utilized  the  fast  Fourier  transform  (FFT)  to  produce 
a slant-range  versus  cross-range  display.  If  the  TTI  system  is  to  be 
effectively  employed  in  an  aerial  warfare  environment  then  real-time 
processing  is  necessary.  In  an  effort  to  speedup  the  signal  processing 
several  alternative  transforms  were  studied  as  possible  substitutes 
for  the  FFT.  The  Karhunen-Loeve,  Cosine  (Sine),  Mellen,  and  Hanhel 
transforms  were  investigated  and  found  to  be  Infeasible  for  use  in  TTI 
Imaging.  T!.e  Walsh  (Hadamard)  transform  was  studied  in  detail  and 
-jested  la  a simulation  program  and  found  that  It  could  not  be  utilized 
In  the  TTI  signal  processing. 

Two  methods  of  converting  from  Che  Walsh  sequency  domain  to  the 
Fourier  frequency  domain  were  studied.  The  first  scheme,  a recursive 
relationship  between  the  arithmetic  and  logical  autocorrelation  functions 
as  presented  by  Robinson  was  dls>'overed  Co  be  incorrect.  The  second,  a 
method  of  computing  the  Fourier  coefficients  from  the  Walsh  coefficients 
of  a function  was  demonstrated  to  be  too  time  consuming  to  be  Implemented 
In  TTI  signal  processing. 

Several  floating-point  FFT  Implementations  were  tested  using  the 
simulation  program.  Also,  several  fixed-point  FFT  algorithms  -were 
derived  and  tested.  All  of  these  were  evaluated  on  the  basis  of  speed  and 
memory  requirements  and  one  fixed-point  FFT  algorithm  was  shown  to  be 
fast  enough  and  accurate  enough  for  Implementation  on  the  TTI  Mini- 
computer . 
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COMPARISON  OF  FAST  FOURIER  TRANSFORMS 
WITH  OTHER  TRANSFOr:MS  IN  SIGNAL 
PROCESSING  FOR  TACTICAL  RADAR 
TARGET  IDENTIFICATION 

I.  Introduction 

The  High  Resolution  Radar  Branch  of  the  Rome  Air  Development 
Center,  Grlfflss  AF£,  New  York  is  interested  in  developing  a 
tactical  target  identification  (TTI)  radar  system  which  can  effectively 
respond  to  the  expected  WARSAW  Pact  threat  against  NATO  in  the  1980-1990 
time  frame.  Radar  target  Identification  of  aircraft  facilitates* 
the  effective  con*‘rol  of  friendly  airborne  Interceptor  and  close  air 
support  aircraft  without  active  Identification  Friend  or  Foe  (IFF) 
devices.  It  also  permits  the  positive  identification  and  determination 
of  uncooperative  or  enemy  aircraft  type  and  mission  without  endangering 
friendly  aircraft  and  allows  the  discrimination  between  actual  enemy 
aircraft  and  decoy  vehicles  (Ref  92). 

The  objective  of  the  TTI  development  program  is  the  generation 
of  two-dimensional  "images"  of  aircraft  by  processing  in  real-time 
the  radar  returns  of  the  aircraft  using  a tactical  pulsed-Doppler 
wideband  radar  system.  Basic  research  is  being  conducted  by  the 
Syracuse  Research  Corporation  (SRC),  Syracuse,  New  York.  The  target 
Identification  simulation  program  used  in  this  thesis  was  written  by  SRC 
and  is  based  on  parameters  modelled  the  same  as  those  of  the  ALCOR*  system 
which  has  been  used  in  actual  radar  Imaging  studies  (Ref  31:1).  In 

*ARPA  (Advance  Research  Projects  Agency) /Lincoln  Laboratory  C-Band 
Observables  Radar. 


essence,  the  two-dimensional  image  can  be  formed  if  a target's  aspect 
angle  changes  sufficiently  relative  to  the  radar  over  some  time 
interval."  The  integration  of  the  received  rader  waveform  with  respect 
to  aspect  angle  will  yield  a slant  range  vs.  Doppler  array  with 
entries  of  radar  cross-section  intensity.  This  array  can  be  displayed 
giving  an  image  of  the  target's  highlights.  The  feasibility  of  such 
an  approach  using  a pulsed-Doppler  wideband  radar  system  and  associated 
signal  processing  has  been  established  by  Rafael  (Ref  68)  and 
Strattan  (Ref  83). 

The  objective  of  this  thesis  is  to  investigate,  evaluate,  and 
-compare  other  orthogonal  transform  (FFT)  currently  being  used  in  the 
signal  processing  portion  of  the  target  identification  system.  The 
goal  is  to  Increase  the  speed  of  the  signal  processing  to  approach 
real-time.  Image  quality,  processing  time,  and  computer  memory 
requirements  must  be  considered  when  Investigating  and  evaluating 
an  alternate  transform. 

The  remainder  of  this  thesis  Is  divided  Into  five  sections. 

Section  II  provides  the  theory,  algorithm,  and  simulation  of  the 
radar  target  identification  system.  Section  III  Investigates  and 
considers  alternate  orthogonal  transforas.  Section  IV  proposes 
two  methods  for  converting  from  the  Walsh  sequency  domain  to  the 
Fourier  frequency  domain.  Section  V compares  several  Fast  Fourier 
Transform  algorithms.  Conclusions  and  recommendations  are  made  in 
Section  VI.  The  Appendices  contain  listings  of  the  computer  programs, 
subprograms,  and  radar  and  target  parameter  data  used  in  this  thesis. 
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II.  Radar  Target  Identification;  Theory, 

Algorithm,  and  Signal  Processing 

This  section  provides  the  background  material  on  Radar  Crosj 
Section  measurement  and  estimation.  The  basic  assumptions  are  given 
and  the  salient  features  of  the  Imaging  algorithm  are  presented. 
Additionally,  the  simulation  program  used  In  this  thesis  Is  explained. 

Target  Radar  Cross  Section 

Radar  cross  section  (RCS),  a.  Is  a measure  of  the  energy  reflected 
from  a target  toward  the  receiving  antenna  (Ref  33:455).  The  RCS  of  a 
target  Is  the  area  assumed  to  Intercept  the  Incident  radiation,  which, 
when  Isotropically  reradlated,  yields  the  actual  power  density  at  the 
receiving  aperture  (Ref  33:455).  This  returned  energy  varies  with 
a multitude  of  parameters  such  as  transmitted  wavelength,  polarization, 
target  geometry,  orientation,  and  reflectivity  (Ref  58:141). 

More  precisely,  the  radar  cross  section  of  an  object  is  proportional 
to  the  far-fleld  ratio  of  reflected  to  incident  power  density,  that  is 

^ _ Power  reflected  back  to  recelver/unlt  solid  angle 
Incident  power  denslty/4ir 

I (Ref  58:142).  For  an  example,  consider  the  RCS  of  a perfectly  conducting 
isotropic  scatterer.  The  power  Intercepted  by  the  radiator  is  the 
'product  of  the  Incident  power  density,  Pj,  and  Its  geometric  projected 
area,  A^.  By  the  definition  of  Isotropic  scattering,  this  power 
is  uniformly  distributed  over  4r  steradians  (Ref  58:142).  For  this 
isotropic  scatterer  then 
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(2) 


Thus,  the  RCS  of  such  an  Tsotioplc  reflector  is  the  geometric  projected 
area  (Ref  58:141). 

For  a .complex  target,  such  as  an  aircraft  or  misslc,  the  RCS  car. 
be  approximated  by  breaking  the  body  into  individual  reflectors 
(scatterers)  and  assuming  that  the  parts  do  not  interact.  In  this  case 
it  can  be  shown  that  the  total  RCS  is  the  vector  sum  of  the  individual 
cross  sections 

M j4ad  2 

a = 1)^  i/ok  exp( — ^— ^)  | (3) 

k=l  ^ 

where  Oj^  is  the  RCS  of  the  k th  scatteror,  dj^  is  the  distance  between 
the  k th  scatterer  and  the  receiver,  and  N the  total  number  of 
scatterers  (Ref  58:144). 

Another  approach  considers  the  relative  phase  angles  between  the 
returns  from  these  N scatterers.  This  approach  leads  to  the  following 
expression  for  the  RCS  of  the  entire  body 

N 2 

CT  = Q'  ^k  exp(j(Ji  ) I 
k=l 

where  ({ij^  is  the  relative  phase  angle  associated  with  the  k th  component 
(Ref  26:974).  It  has  been  shown  that  the  RCS  of  a target  is  related 
to  the  frequency  response  of  the  object,  Gfjw) , by  the  following 
relationship, 

o = icao))!^  (5) 

o 

(Refs  51:1651;  83:5).  Tne  RCS,  a , can  be  thought  of  as  the  spectrum 

s 

of  the  complex  target.  In  general,  G(ja))  will  be  aspect  angle  dependent 
except  for  a spherically  symmetric  object  (Ref  51:1651). 

Now,  if  a CW  or  pulses  radar  signal  is  reflected  by  a target  moving 


4 


at  a velocity,  v^,  relative  to  the  radar  receiver,  the  whole  spectrum 
would  be  translated  In  frequency  by  the  Doppler  shift,  fjj,  where 


(6) 


where  c Is  the  propagation  velocity  (speed  of  light) , f^  Is  the  trans- 
mitted carrier  frequency,  R or  v^  Is  the  range  rate  or  radial  velocity, 
and  fjj  is  the  Doppler  8l:‘'t  (Refs  58:5;  47:357). 

In  a wideband  pulsed-Doppler  radar  system,  the  received  Doppler 
frequency  spectrum  Is  considered  the  target's  radar  cross  sectlop  which 
Is  aspect  angle  dependent  (Ref  58:173).  Figure  1 shows  the  effect  of 
return  power  from  an  aircraft  where  the  surface  scatterers  making 
up  the  composite  target  echo  can  create  a transition  from  phase  addition 
to  phase  cancellation  and  change  the  cross  section  drastically 
(Ref  33:26).  The  spectra  of  RCS  flucuatlons  can  be  described  In  terms 
of  several  effects  with  the  airframe  the  most  Important  contrlbuter 
(Ref  58:173).  The  airframe  spectrum  Is  due  to  the  relative  motion 
between  the  various  scattering  points  on  the  fuselage  and  wings.  This 
relative  motion  occurs  as  the  aircraft  aspect  changes  (Refs  58:173; 
51:1651;  26:973). 

llie  resulting  spectral  width  Is  proportional  to  the  transmitted 
frequency  (Refs  58:73;  83:3).  The  frequency  domain  will  give  a band- 
width, therefore,  of 

B • 1/T  (7) 


where  T is  the  pulse  duration  length  (Ref  83:4) 


Figure  1.  Return  Power  from  an  Aircraft  as  a Function 
of  Azimuth  Angle 


RCS  Measurement 

The  measurement  of  the  radar  cross  section  of  a target  makes  use 
of  Fourier  transform  theory  and  it  will  be  shown  that  |G(ju)p  Is  In 
reality  the  Power  Spectral  Density  (PSD)  or  Power  Spectrum  of  the 
Impulse  response,  g(t),  of  the  target. 

Fourier  Integral.  The  Fourier  Integral , defining  a Fourier 
transform  pair.  Is  given  for  f(t)  specified  on  the  Interval  (-»,»)  as 

F(ju)  - f (t)exp(-j(ot)dt  (8) 
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snd 


F(ju)exp(J&)t)du 


Power  Spectral  Density.  If  f(t)  Is  specified  on  the  Interval 
(-*••)  and  If  f(t)  and  F(j(ii)  are  a Fourier  transform  pair,  then,  the 
Power  Spectral  Density  of  the  function,  f(t)  Is  defined  as  the  abso- 
lute value  squared  of  the  Fourier  transform  of  f(t).  If  P(u)  Is  the 
Power  Spectrum,  then 


P(u)  - |F[f(t)]|2 
- |F(Jw)P 


where  f[*]  Is  the  Fourier  transform  operator.  Also,  If  the  autocor- 
relation function  of  f(t),  R(t),  Is  Fourier  transformed,  then  the 
following  relationships  will  be  found  to  be 


f 

R(t)  - y J f (t)f(t+T)dt 
-T/2 


- f(t)  * f(-t) 


where  * denotes  convolution,  but 


F[f(t)*f(t)]  - F(Jo))*F(-ju) 
- |F(jw)|2 


\ 

\ / 


80  that 

p(u)  - F[R(t)1  (13) 

Equation  (10)  is  known  as  the  direct  cethod  for  obtaining  the  power 
spectrum  of  a tine  series  and  Is  equal  to  the  Fourier  transform 
squared  of  the  function  (Ref  76:14,15).  Equation  (13)  Is  known  as 
the  direct  method  or  the  Wlener-Khlntchlne  theorem.  It  states  that 
the  power  spectrum  Is  the  Fourier  transform  of  the  autocorrelation 
function  (Refs  76:3;  28:128). 

Therefore,  It  can  be  seen  that  Equation  5 Is  In  fact  the  Power 
Spectral  Density  of  g(t),  the  Impulse  response  of  the  target.  (That 
which  would  be  measured  at  the  Doppler  filter.) 

In  this  section  the  power  spectrum  Is  defined  for  an  Infinitely 

long,  continuous  time  function,  £(t)  or  g(t).  However,  In  practical 

) 

situations,  only  a finite  amount  of  time  Is  even  available  to 
observe  the  time  function,  particularly  If  a monostatic  radar  system 
Is  used. 

Since  the  signal  Is  In  effect  truncated,  the  effects  on  che 
Power  Spectral  Density  resulting  from  the  truncation  of  the  data  set 
must  be  considered. 

RCS  Estimation 

Since  the  RCS  of  a target  cannot  be  completely  determined.  It 
must  be  estimated.  The  Discrete  Fourier  transform  (DPI)  can  be 
used  to  compute  an  estimate  of  the  Radar  Cross  Section  of  a complex 

target. 

Discrete  Fourl  r Transform.  Let  f(n)  be  defined  by  N samples. 

The  Discrete  Fourier  transform  pair  are  defined  as 
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N-1 

F(k)  = 1 f (n)exp(-j27r(kn)/N)  (14) 

n=0 


and 


N-J. 

f(n)  = F(k)exp(j27r(kn)/N)  (15) 

k=0 

where  the  exponential  function  is  periodic  of  period  N (Ref  60:100), 
Discrete  Power  Spectral  Density.  The  sample  power  spectral  den- 
sity function  or  what  is  known  as  the  raw  perlodogram  Is  the  DFT  of 
R(t),  the  autocorrelation  function  of  f(n)  (Ref  88:13).  Define 

1 N 

p(oj)  ^ R(T)exp(-jTu)  (16) 

t=-N 


where  T Is  an  inteser. 

Or  conversely,  the  perlodogram  can  be  calculated  as  the  modulus  of 
the  DFT  of  f (n) . 


N-1 

P(a))  - ■ I f (n)exp(-jmo/N) 

I »^2TrN  n=l 


(17) 


Equation  (16)  is  analogous  to  the  continuous  Wiener-Khintchlne  theorem 
only  in  discrete  form. 

As  it  turns  out,  whichever  way  is  used,  the  raw  perlodogram  is 
an  unsatisfactory  estimate  of  the  power  spectrum  unless  the  signal 
is  perfectly  periodic  and  noiseless.  Therefore,  the  perlodogram  of 
a stochastic  process  will  be  an  unstable  estimate,  erratic  In 
appearance  and  behavior  (Ref  88:14).  The  variance  of  the  fluctuation 
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pw 


of  the  perlodogram  about  the  true  power  spectrum  does  not  decrease 
to  zero  as  M approaches  infinity,  as  it  should  for  a well  behaved 
estimator  (Refs  88:14;  30:109).  Its  variance  is  independent  of  N, 
and  the  probability  distribution  of  the  sample  perlodogram  is  a Chi- 
squared  distrlbutinn  with  two  degrees  of  freedom  (Ref  88:80).  To 
bring  the  variance  down  and  stabilize  the  estimate,  it  is  necessary 
to  do  some  form  of  spectral  averaging  (Refs  76:13;  88:82;  62:548). 

Currently,  the  "practice"  of  power  spectral  estimation  has  a 
strong  empirical  basis  because  most  optimum  techniques,  such  as 
timTcimiim  likelihood  esClnuitlon,  require  more  information  about  the 
signal  t^an  is  usually  available  (Rev  62:532).  As  a result  trade-offs 
are  involved  between  different  techniques  such  that  there  is  no 
general  agreement  on  the  best  method.  The  reader  is  directed  to 
the  literature  for  a more  detailed  look  at  power  spectral  estimation. 
Davenport  and  Root  (Ref  28),  Welch  (Ref  90),  Webb  (Ref  88),  Sentman 
(Ref  76)  and  others  are  excellent  references.  Oppenhelm  and  Shafer 
(Ref  62)  provide  an  excellent  overview  of  the  power  spectral  estima- 
tion problem.  Deutsch  (Ref  30)  Is  a recommended  reference  for 
estimation  theory.  Blacksmith,  et  al,  (Ref  17)  contains  an  extensive 
bibliography  on  work  done  on  radar  cross  section  measurement. 

Estimating  Discrete  Power  Spectrums 

Let  f(t)  be  defined  on  the  interval  (-“,»).  If  f(t)  is  truncated 
by  multiplying  by  a data  or  observation  window,  w(t),  (t  can  be 
considered  either  continuous  or  discrete)  such  that 


10 


w(t)  - < 


t!  < T/2 
tj  > T/2 


(18) 


1 

0 


then  the  truncated  data  set  becomes 


hu)  ” f(t)*w(t) 


(19) 


as  shown  In  Figure  2 on  the  next  page. 

The  functlcn  h(t)  now  represents  the  truncated  data  set  available 
from  which  the  power  spectrum  Is  to  be  calculated.  The  Power 
Spectral  Density,  P__(m)  of  h(t),  representing  the  apparent  power 

ap 

spectrum  of  f(t},  Is  then 


Pap(“)  - |F[f(t)*w(t)]|2 
■ |F(Ju)*W(ju)|2 
- lF(jaj)|2*  |W(ja))|2 


(20) 

(21) 

(22) 


Thus,  the  PSD,  |F(jb))|2,  Is  modified  by  a convolut  ion  with  |w(ja>)|2, 
the  Fourier  transform  of  the  data  window.  K(ju))  Is  called  the 
frequency  window  and  Is  a sin  x/x  (sine  x)  function  as  shown  In 
Figure  3. 

Convolution  by  the  frequency  window  causes  a certain  degree  of 
smoothing  In  the  calculated  PSD,  but  a small  amount  of  leakage  via  the 
sldelobes  from  nearby  frequency  bands  Into  the  frequency  band  of 


11 


figure  2.  Truncated  Data  Set 


interest  (Ref  76:12).  Rormally,  this  Is  not  serious,  but  If  there 
are  large  veil  defined  peaks  In  the  power  spectrum  and  when  convolved 
with  the  frequency  ' indow,  they  act  to  reproduce  the  window,  producing 
spurious  peaks  In  the  PSD  corresponding  .:o  the  sldelobes  of  the 
frequency  window.  These  spurious  peaks  may  be  mistakenly  identified 
as  structural  details  of  the  true  power  spectrum  when  In  fact  they 
are  merely  artifacts  created  by  the  truncation  of  the  original 
data  set  (Ref  76:13). 

Several  methods  exist  for  sldelobe  suppression  Three  common 
methods  Include  data  tapering,  sectioning  and  averaging,  and  data 
weighting  (Refs  88:84;  2:548;  90:56).  The  simplest  and  most 

straight  forward  method  Is  the  use  of  data  weighting.  Several  types 
exist  and  the  most  widely  used  Is  the  Hamming  weight  function 


L 


2 


(Sine  x) 


Figure  3.  Sin  x/x  (Sine  x)  Function 

I 

I 

i 

WgCn)  < 0.54  - 0.46co8(2iin/(N-l)).  (23) 

for 

0 < n < N-1 

where  N Is  the  total  number  of  samples  (Ref  (62:241-242). 

The  statistical  reliability  of  PSD  estimates  Is  discussed  by 
Webb  (Ref  88),  Sentman  (Ref  76),  and  Davenport  and  Root  (Ref  28). 

Assumptions 

The  basis  of  the  Imaging  technique  used  In  this  thesis  Is  based 
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on  the  assumption  that  the  body  being  observed  is  rigid,  and  coasists 
of  sections  that  behave  as  point  scatterers  (Ref  56:1-4).  It  is  also 
assumed  that,  relative  to  the  observer,  the  object  has  some  motion  about 
a fixed  set  of  axes  (Ref  56:1-4).  Figure  4.  shows  a typical  body  with 
the  axes  centered  at  the  center  of  rotation  of  the  body  and  with  the 
rotation  vector  normal  to  the  x-y  plane  (Refs  56:1-4,4;  68:5).  The 

body  l3  assumed  to  remain  in  one  range  bin.  'Hie  separation  between 
transmitted  pulses  is  sufficient  to  prevent  overlapping,  which  is  met 
by  the  following  condition 

4T,  - N T.  (24) 

d r h 

where  T.  is  the  duration  of  the  total  waveform,  T,  is  pulse  duration 
within  the  waveform,  and  is  codelength  required  for  PRF  staggering 
(Ref  82:4).  The  maximum  unambiguous  range  of  the  radar  is  determined 
by  the  burst  length,  T^,  so  that 

■ '’•d/2 

and  the  first  Doppler  ambiguity  is  sufficiently  removed  to  correspond 
to  double  velocity  of  the  fastest  target  to  eliminate  foldover,  or 

/-=  C26) 

where  T^  is  the  average  Interval  between  subpulses  (Ref  83:4). 

Additionally,  it  is  assumed  the  signe.'  processing  is  of  first-order, 
which  means  that  integration  Intervals  are  short  enough  for  certain 
ifi^earltles  to  pervall,  that  is,  the  body  has  linear  motion  (no  acceler- 
ation) during  the  processing  Interval  in  which  the  target  is  rotating 
' ■ It  some  effective  rotation  center  (Ref  68:3).  Iii  other  words, 

thr.  Doppler  shift  produced  by  the  rotating  scattering  points  la 
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Arbitrary  location  but  point  rotater  about 
fixed  x-7  center  with  radius  L 

Two-way  Doppler  shift  f^  : V^/(A/2) 

■ Radial  Velocity 

V ■ L ■ Total  velocity  of  point 

n ■ rotation  rate 
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V ■ V sin  0 ■ n L sin  0 
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L sin  0 ■ X ■ position  along  X“axls 
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X - f .(X/2) 
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Figure  4.  Range-Doppler  Image  Relationships 
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considered  approximately  constant  during  a signal  proces'slng  interval. 
The  rotation  can  be  viewed  as  a changing  aspect  between  the  radar  line 
of  sight  (LOS)  and  the  target. 

Imaging  Algorithm 

The  radar  imaging  simulation  program  models  the  target  as  a 
collection  of  scattering  points  p^,  with  each  point  assigned  an  RCS 
v’hlch  is  assumed  to  be  aspect  angle  Independent  (Befs  31:1;  57:2). 

At  equal  time  Increments  the  RCS  is  modelled  as  a function  of  frequency 
according  to 


N ^ ^ 

I *^a~exp  C-jAirf  k'f  |c) 
-1  ^ 


(27: 


where  N represents  the  number  of  scatterers,  k is  the  RLOS  direction 
vector,  and  r^  is  the  range  from  an  origin  on  the  body  to  the  point 
p^  (Ref  31:1).  For  one  image,  "flight"  continues  until  the  body 
has  undergone  sufficient  aspect  angle  change  to  give  a cross-range 
resolution  of  1.5  feet  or  approximately  2.9  degrees  of  change,  as 
determined  by 


AX 


fDCX/2) 


an) 

A6m 


(28) 


where  AOm  is  the  aspect  angle  change  of  the  body  relative  to  the  RLOS  and 
is  the  effective  vehicle  rotation  rate  (angular  velocity)  (Ref  68:5). 
The  range  resolution  is  shown  in  Figure  5 is  determined  by  the 
radar  bandwidth.  The  dwell  time  on  target  needed  to  produce  an  image  with 
a given  AX^  is 

Wii  ■ '>•“  O’) 
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where  it  can  be  shown  that 

sin  0 (30) 

where  (2^  is. the  Instantaneous  RLSO  angular  velocity  due  to  the  varying 
aspect  between  the  target  track  and  the  RLOS  and  sin  6 Is  the 
Instantaneous  velocity  of  the  a ^le  perpendicular  to  the  RLOS 
(Refs  68:6;  56:1-5). 

Since  a moving  target  has  many  degrees  of  freedom  of  motion,  the 
effective  vehicular  rotation  may  be. brought  about  by  changing  aspect 
between  the  RLOS  and  the  target  track,  and  by  changes  In  aspect  due 
to  target  motion  abouts  Its  center  of  mass  (Ref  68; 6). 

The  magnitude  of  0^  Is  used  to  scale  the  image  from  Doppler  to 

cross-range  units,  and  the  Image  projection  plane  Is  the  plane 

-»•  

perpendicular  to  0^(Ref  68:6).  This  assumes  that  0^  Is  essentially 

constant  In  both  magnitude  and  direction  over  the  processing  Interval 

(Ref  68:6). 

Signal  Processing 

Data  from  a total  time  span  T,  where  total  angular  change  Is 
given  by  ^ 

A0J,  - V 1 (31) 

I 

and  a range  window  R centered  on  the  target  are  sampled  and  collected 
and  stored  In  an  array.  Both  the  amplitude  and  phase  of  the 
are  stored.  The  sampling  Increment  In  range  Is  R and  In  time 
where 

At  « 1 /radar  pulse  repetition  frequency 


(32) 


To  put  the  data  Into  the  assumed  form  of  a collection  of  points 
rotating  about  a rotation  center,  it  is  necessary  to  compensate 
for  the  motion  of  the  rotation  of  the  rotation  center  with  respect  to 
the  radar  since  if  the  target  is  rotated  slightly,  the  phase  of  the 
center  will  change  (Ref  68:7).  This  can  be  done  by  obtaining  the 
distance  to  the  rotation  center  as  a function  of  time  with  radar 
tracking  data,  and  then  by  subtracting  the  calculated  signal  phase 
at  each  time  from  all  of  the  phase  values  in  a pulse  return  recorded 
at  that  time.  Alternately,  the  phase  recorded  for  an  actual  discrete 
target  point  in  the  signal  may  be  used  as  a reference  for  compensating 
all  the  phases  in  the  return  (Refs  56:1-7;  68:7).  This  process  is 

known  as  aligning  or  "cohering"  the  data.  Next,  the  data  are  Hamming 
weighted,  discussed  earlier,  for  sidelobe  control  and  Fourier 
transformed  along  constant  slant  range  lines  to  produce  the  image  output 


(Ref  68:7). 

The  unambiguous  cross-range  window  Is  given  by 

Y . 

*ACR  A0j^ 

The  corresponding  unambiguous  range  window,  the  Doppler 


(33) 


window,  are  simply 


-Mi 

ACR  At 


ft/s 


D 


A(K 


with  X in  feet  and  At  in  seconds  (Ref  68:8), 

The  cross-range  grid  Increment,  AX  is  given  by 

AX  = X.„_/Cnumber  of  output  points  in  the 
transform) 


(34) 

(35) 


(36) 
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In  general,  AX  < AX^,  the  cross-range  resolution,  in  order  that  the 
sampling  not  be  too  coarse  to  miss  significant  output  features  (Ref  68:9). 
Strattan  (Ref  83)  points  out  that  the  number  of  profiles  (HCS  estimates) 
Integrated  should  be  at  least  equal  to  the  ratio  of  the  maximum  cross-^ 
range  dimension  of  the  target  to  AX^  and  that  the  transforms  should  be 
performed  at  slant-range  Intervals  no  larger  than  the  range  resolution 

AR  » c/2B  (37) 

where  B * 1/T,  the  Doppler  bandwidth  of  the  system.  The  angular  Interval, 
needed  for  the  cross-range  scale  factor  may  be  determined  approximately 
from  flight  path  tracking  data.  Figure  5 shows  the  resolution  relation- 
ship necessary  for  good  Imaging  (Ref  68:5). 

Detailed  mathemetlcal  first-order  range-Doppler  processing  is 
given  by  Rafael  (Ref  68). 

Simulation  Program  TGTID 

The  Simulation  Program  TGTID  (Appendix  A)  used  In  this  thesis  was 
written  by  researchers  at  the  Syracuse  Research  Corporation, 

Syracuse,  New  York.  The  listings  of  the  program  and  Its  supporting 
subroutines  are  located  In  the  Appendices. 

The  "input  data"  Is  coherent  radar  cross-section  data  as  a 
function  of  frequency.  It  Is  Inversed  Fourier  transformed  to  radar 
cross-section  as  a function  of  range.  This  data  Is  first  aligned 
(cohered)  so  that  the  first  peak  of  each  range  sample  occupies  the  same 
range  bln  and  then  a phase  adjustment  Is  made  giving  the  first  peak 
zero  phase. 

For  fixed  range,  the  adjusted  data  (now  radar  cross-section  as  a 
function  of  pulse  repetition  Interval)  Is  Fourier  transformed  along 
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Figure  5.  Range-Doppler  Resoletlon  Relationships 


constant  range  cells  to  give  the  RCS  In  terms  of  Doppler  frequency 
which  Is  then  scaled  to  cross-range.  The  result  is  a cross-range 
versus  range  "radar  Image"  \diose  entries  are  the  associated  RCS*s 
(Ref  56:2).  This  array  Is  then  displayed  and  "squared  up"  so  that 
an  undlstorted  "picture"  of  the  aircraft  RCS  response  is  realized.  A 
block  diagram  of  the  processing  Is  shown  in  Figure  6 (Ref  68:8). 
Program  TGTID  Description 

The  main  program  TGTID  (Appendix  A)  reads  in  the  radar  (location 
and  frequency)  and  target  (scatterer)  parameters.  Subroutine  CINIT 
(Appendix  B)  computes  the  number  of  samples  and  order  of  these  samples 
as  a power  of  two  based  on  the  given  parameters.  Subroutines  SLARTV* 


(Appendix  B),DOTP  (Appendix  B) , CTRAN  (Appendix  B),  FFT  (Appendix  C) , 


and  ROLL  (Appendix  B)  are  used  to  "generate"  the  RCS  data  as  a function 
of  range.  Subroutines  CIMAGE,  (Appendix  B) , HAMI^GT  (Appendix  B), 

FFT,  and  ROLL  are  used  to  process  the  data  and  Subroutines  FLOTID 
and  BUFOUT  (Appendix  G)  are  used  to  display  the  synthetic  Image 
created  from  the  RCS  data.  A complete  description  of  each  subroutine  Is 
given  with  their  listing  In  the  Appendices.  Radar  and  Scatterer 
parameters  used  In  this  thesis  are  In  Appendix  H. 

For  this  thesis,  four  scatterers  were  simulated  and  their  returns 
processed.  A representation  of  the  target  (four  scatterers)  and 
the  radar  set  Is  shown  In  Figure  7.  A sample  "Image"  of  the  four 
scatterers  as  "processed"  by  the  Syracuse  Research  Corporation's 
computing  system  using  this  simulation  program  with  FFT6  Is  shown  In 
Figure  8 and  was  used  as  a standard  of  comparison  for  the  alternative 
forms  of  processing  used  In  this  thesis. 


Z 


Figure  7.  Simulated  Target  and  Radar  Configuration 
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Subroutine  PLOTID  (Appendix  G)  was  written  to  display  the 
generated  "image"  uslnp  a CALCOMP  plotter.  Each  display  is 
generated  by  comparing  each  element  in  the  image  matrix  and  plotting 
a symbol  if  the  value  exceeds  a threshold  level.  In  this  way,  only 
the  prominent  power  spectral  pealer  or  highlights  of  the  scatterers 
are  plotted.  Figures  9 through  14  show  the  images  plotted  using  the 
Simulation  Program  TGTID  with  Subroutine  PTOTID  and  Subroutine  FFT6. 
The  threshold  set  in  Figures  9 through  11  are  set  too  "low".  The 
threshold  set  in  Figures  12  and  13  are  in  the  proper  range,  where  the 
threshold  set  in  Figure  14  is  too  "high".  The  threshold  level  in 
an  actual  system  would  be  set  dynamically  with  range  information  and 
signal  "strength". 

All  programming  is  done  in  FORTRAN  IV  Extended  on  the  CDC  6600 
CYBER  74  System.  The  plots  are  generated  using  a CALCOMP  Plotter. 


III.  Evaluation  of  Orthogonal  Transforms 


The  use  of  orthogonal  transforms  is  Investigated  in  this  Section. 

A transform  usually  possesses  some  attributes  that  make  it  desirable 
to  use  in  a particular  application.  These  characteristics  Include; 
applicability  to  the  problem  at  hand;  mean  square  error,  probability 
of  error,  computational  advantage,  programmable  in  some  computer 
language,  and  computer  memory  required  to  compute  the  transform. 

These  characteristics  are  often  used  as  criteria  for  determining  the 
acceptability  of  the  transform,  especially  for  implementation  on  a 
digital  computer. 

The  Karhunen-Loeve , Cosine  (Sine),  Kellln,  and  Hankel  transforms 
‘are  Investigated  and  found  not  to  be  suitable  for  this  application. 

The  Walsh  (Hadamard)  transform  is  examined  in  detail  because  of  its 
slfflllaratles  vlth  the  Fourier  transform.  Its  computational  ease 
and  speed,  and  the  existence  of  an  analogous  Wlener-Khlntchlne  theorem. 
Karhunen-Loeve  Transform 

The  optimum  transform  for  data  compression  and  for  satisfying  the 
minimum  mean  square  error  criteria  i?  the  Karhunen-Loeve  transform 
(Ref  9:123).  Its  most  common  application  is  found  in  image  and  picture 
transform  coding  where  data  compression  is  highly  desired  (Ref  48:64,  65) 
The  transform  is  composed  of  eigenvectors  of  the  correlation  matrix  of 
the  original  signal,  picture,  or  class  of  Images  to  be  coded  (Ref  9:124). 
The  reader  is  referred  to  Andrews  (Ref  9)  for  a detailed  treatment  of 
the  Karhunen-Loeve  transform. 

There  are  two  major  problems  associated  with  the  use  of  the 
Karhunen-Loeve  transform.  The  first  is  that  a great  amotut  of 
computation  must  be  performed  (Refs  9:125;  10:41).  The  correlation 
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matrix  must  be  estimated  If  It  Is  not  known.  Next,  the  correlation 

matrix  must  be  diagonalized  to  determine  Its  eigenvalues  and  eigenvectors. 

Finally,  the  transform  Itself  must  be  taken.  In  general,  there  la  no 

fast  computational  algorithm  for  the  transform.  Since  It  Is  usually 

not  separable  (Ref  9j125)  fFor  example:  for  N data  points,  the 

2 2 

computatlona}  load  Is  2(n  ) ; for  128  points,  tnls  Is  over  500 
million  multiplications,  which  Is  too  many  to  be  feasible.)  The 
second  difficulty  Is  that  the  mean  square  error  Is  not  a valid  error 
criterion  for  most  applications  (Ref  9:125)  Including  this  one. 

Cosine  (Sine)  Transform 

The  Cosine  (Sine)  transform  has  been  found  to  compare  favorably 
with  that  of  the  optimal  Karhunen-Loeve  transform  (Refs  4:90;  48:71,72). 

The  Discrete  Cosine  Transform  (DCT)  of  a data  sequence  xCn), 
n “ 0,1,2, ..., (n-l)  Is  defined  as 

G,(0)  - I xtn)  (38) 

* " n-0 

G Ck)  - I f ^ x(xi)co8  (39) 

* " n^O  2N 

k « 1,2,..., (N-l) 

where  N Is  the  total  number  of  samples  and  G^(k)  Is  the  k th  DCT 
coefficient  (Ref  4:90).  Ahmed  et  al.  state  that  Equation  39  can  be 
expressed  as 

, 2N-1  . 

G/k)  - jRe  {exp((-jkn)/2N)  J x(n)Wjj‘“}  UO) 

n»0 

where  W exp(-j2n/2N) , j - i^,  x(n)  ■ 0 for  m - N,  (N+1),..., 

(2N-1),  and  Re{*}  Implies  the  real  part  of  the  term  enclosed  (Ref  4:91). 
Fton  Equation  40  it  follows  that  the  N DCT  coefficients  can  be  computed 
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using  a 2N-polnt  fast  Fourier  transform  (^ef  4;91).  Similarly,  If  a 
discrete  Sine  transform  were  desired,  the  Se{*}  In  Equation  40  would 
be  replaced  by  Im{ * } , which  denotes  the  Imaginary  part  of  the  term 
enclosed.  It  can  be  seen  that  the  computational  speed  of  the  DCT  or  DST 
Is  slower  as  that  of  the  FFT  since  twice  as  many  points  must  be 
transformed.  No  other  fast  transform  algorithm  curren'ly  exists  for 
the  DCT  or  the  DST. 

Mellln  Transform 

The  Mellln  transform  possesses  the  unique  property  of  "scale" 
Invariance  (Rsf  22:78).  That  Is,  scale  changes  In  the  Input  do  not 
produce  scale  changes  In  the  output.  The  Melllir  transform  M(u)  of  a 
function  f(x)  la  defined  by 

M(u)  - I f(x)  x’“"^dx  (40) 

0 


and  the  Inverse  transform  Is  given  by 
£(x) 


_JL  r"" 

" Jy-la 


M(u)  X '*du 


(41) 


where  y Is  chosen  so  that  the  Integral  exists  (Ref  83:13).  The  discrete 
Mellln  transform  Is  given  by 


M(kAu)  - I fClAx)ClAx)^’^“"^  Ax  (42) 

1-1 

where  N la  the  number  of  samples  of  f(x),  and  the  Input  and  transform 
space  resolutions  are  Ax  and  Au  respectively  (Ref  22:79). 

No  fast  computational  algorithm  has  yet  been  found  for  Implementing 
the  Mellln  transform  directly.  Processing  time  for  a digital  computer 
can  be  quite  long  (Ref  22:80).  Cesasent  and  Psaltls  point  out  that  a 
digital  Mellln  transform  can  also  be  realized  by  exponantally  sampling 
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the  input  and  then  performing  an  FFT  on  this  data  (Ref  22:78,80),  It 
is  concluded,  therefore,  that  the  Mellln  transform  offers  no 
advantage  over  the  FFT  processing  currently  being  performed  in 
the  simulation  program.  For  a further  comment  on  this  subject  -see 
Chapter  VI. 

Hankel  Transform 

The  Eaakel  transform  is  useful  if  symmetry  exists  about  an  axis 
and  if  polar  coordinates  are  appropriate  (Ref  85:46).  The  Hankel 
transform  pair  are  defined  by 


F(10  - f(x)  J Ckx)  dx 


£(x)  - FCk)  J (xk)  dk 


vhere  J (kx)  is  the  Bessel  function  of  the  first  kind  of  order  n 
n 

(Ref  85:46).  This  application  can  be  seen  to  demonstrate  some  elements 
of  symmetry  about  an  axis,  say  the  aspect  angle  of  the  body  at  T/2  of  the 
observation  Interval,  but  no  method  to  convert  the  received  Doppler 
waveform  to  polar  coordlnnces  could  be  found.  Additionally,  this 
transform  requires  the  calculation  of  certain  boundary  conditions 
which  is  a very  difficult  problem  (Ref  85:46).  It  was  concluded  that 
the  Hankel  transform  could  not  be  applied  to  this  problem. 

Walsh  Functions  and  Walsh  Transforms 

Walsh  functions  are  a complete  orthonormal  set  of  square  wave 
functions  that  are  finding  increasing  use  in  various  digital  signal 
processing  applications  (Ref  29:137,  Ref  44).  They  exhibit  similarities 
to  the  trigonometric  sine  and  cosine  functions  in  many  of  their 
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properties  (Refs  87;  4;  5).  The  bivalued  characteristic,  orthogonality, 
and  cciTiputational  advantages  of  these  functions  are  the  basis  and 
motivation-  for  'iheir  detailed  study  in  this  section. 

Definition.  Continuous  Walsh  functions  may  be  defined  in  several  ways 
(F.af  52:211).  They  may  easily  be  defined  as  products  of  Rademacher 
finctions  (Refs  j2:212;  84:4). 


The  Rademacher  functions  (Refs  29:177;  84:4)  are  defined  by 


Ro^e) 


1.  0 1 0 

-1,  1 < 0 £ 1 

2 


Rq(0  + 1)  - RqC0) 


Rj^CQ)  - RqC2"0),  n - 1,2,3,...  . 


Figure  15  shows  the  first  five  Rademacher  functions  (Ref  55:39). 

To  form  the  Walsh  function  wal(n,0),  first,  form  the  binary 
representation  of  n,  then  form  the  Gray  code  version  of  n,  and  multiply 
together  Rademacher  functions  according  to  the  1 bits  in  the  Gray  code 
(Ref  52:212). 


If  n in  binary  is 


N = b b -b  , ...  b 
m m-1  m-2  o 


then  n in  Gray  code  is 


Xvhsre 


^ 9 • t * 8 

m m— i m— i o 


m ra 
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where  9 is  modulo-2  addition  with  no  carries  (Ref  52:212).  The  Walsh 
functions  can  now  be  defined  by 

walCO.e)  ^1  ^22) 

wal  (n,e)  - R ®“(0)  R . . . R *°(e)  (53) 

m m-i  o 

where  n is  the  order  and  6 is  normalized  time  (Ref  84:4).  This  can  also 
be  written  as 

n 

wal(n,e)  - I g/KT.(e)  (54) 

k-0  “ 

where  the  summation  symbol  denotes  modulo-2  summing  (Ref  29:187). 

For  example: 


- IIIO2  - 101 

(55) 

wal(6,e}  - R2(0)R^°(0)RQ^Ce) 

(56) 

- R2C0)RoC0) 

(57) 

This  procedure  to  find  a particular  Walsh  function  is  easily 
remembered.  It  should  be  noted  that  the  first  argument  of  a Walsh 
function  denotes  its  "sequency”.  Sequency,  as  defined  by  Harmuth, 
is  the  number  of  zero  crossings  ur  sign  changes  of  the  Walsh  function 
in  the  half-open  interval  (0,1)  (Ref  44:50). 

Harmuth  uses  the  notation  Wal(j,G)  to  define  the  Walsh  function  and 
further  defines 


37 


WalCQ,0)  - 1 

(58) 

Cal(n,e)  - Wal(2n,6) 

(59) 

Sal(n,0)  - Wal(2n-l,e) 

(60) 

where  the  Cal  and  Sal  namea  are  used  to  parallel  the  cosine  and  sine 
functions  (Rei.  44:22).  Several  notation  schemes  are  used  In  the 
literature  and  are  summarized  by  Meek  (Ref  55*11).  The  first  eight 
continuous  Walsh  functions  are  shown  In  Figure  16. 

Discrete  Walsh  Functions.  Discrete  Walsh  functions  are  sampled 
versions  of  the  continuous  set  (Ref  77:457).  Shanks  assumes  that  the 
discrete  functions  are  Infinite  In  extent,  and  are  periodic  with  period 
N,  where  N Is  an  Integral  power  of  two  (Ref  77:457).  Thus  a complete 
orthogonal  set  will  have  N distinct  functions,  designated  as  wal(n,m). 

The  complete  set  Is  represented  over  the  range  n ■ 0,1,...,N-1  and 
m <"  0,1,...,N-1.  The  first  two  discrete  Walsh  functions  are  defined 
aa 

wal(0,m)  “ 1,  n ■ 0,1,2,. ..,N-1 

wal(l,m)  ■>  l,m  ••  0,1,2, ...,  (N/2)-l 

- -l,m  - N/2,  (N/2)  + 1,...,N-1  ’ 

Various  Iterative  equations  have  been  used  to  generate  the  remainder  of 
the  set,  but  Henderson's  seems  to  be  the  most  convenient 

wal(n,m)  - wal([N|2],2m)  . wal(N-2 [N/2] ,m)  (64) 

where  [N/2]  Indicates  the  Integer  part  of  N/2  (Refs  45:51;  77:457). 


(61) 

(62) 

(63) 
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Shanks  shows  that  this  recursive  equation  generates  a complete 
orthogonal  set  (Ref  77:459).  Figure  17,  shows  the  first  eight  discrete 
Walsh  functions  of  length  eight  generated  by  Equations  61,  62,  and  63 
(Ref  77:457). 

The  Walsh  functions  may  also  be  represented  In  matrix  notations. 

Let  N 2 , where  k Is  a positive  Integer,  then  the  N th  order  Walsh 
matrix  Is  constructed  by  sampliig  the  first  N Walsh  functions  once 
In  N equal  subintervals  of  (0,1)  (Ref  87:5).  The  matrix  constructed 
from  sampling  the  sequency  ordered  Walsh  functions  (Figure  16)  results 
In  a sequency  ordered  NxN  Walsh  matrix.  Figure  18a,  shows  the  sequency 
ordered  Walsh  matrix  for  N ••  8. 


' V 


A second,  simpler,  method  which  generates  Walsh  matrices  In  what  Is 
termed  "natural"  order  Is  a Kroneker  product  of  p second  order  Walsh 
matrices  (Refs  85:5;  78:177).  This  form  or  the  Walsh  matrix  Is  of::en 

referred  to-  as  the  "Hadamard"  matrix.  A second  order  Walsh  matrix  Is 
defined  by 


W2-H2 


#* 

1 1 
,1  -1 


(65) 


Hadamard  matrices  of  higher  order,  for  N a power  of  tv;o,  are  generated  by 
the  Krontker  product  operation,  such  that 


■«N 

Hh  - 


(66) 


Figure  18b,  shows  the  natural  ordered  Walsh  matrix  or  Hadamard  matrix  for 
n ■ 8. 

Both  the  Walsh  matrix  and  the  Hadamard  matrix  are  square  arrays, 
whose  rows  and  columns  are  orthogonal  to  each  other,  that  Is,  the 
product  of  the  matrix  and  Its  transpose  Is  the  Identity  matrix  times  N, 
where  N Is  the  order  of  the  matrix  (Ref  78:177). 

H'h’^  - N*I  (67) 

Walsh  Transform.  Since  the  Walsh  functions  form  a complete, 
orthonormal  set  over  the  interval  (0,1),  any  absolutely  integrable  function 
defined  over  the  Interval  can  be  expanded  into  a series  of  Walsh 
functions  analogous  to  the  Fourier  expansion  of  such  a function  (Ref  44:45). 
The  discrete  Walsh  transform  of  a function  is  defined  by 
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and  the  Inverse  transform  Is  given  by 

N-1 

f(n)  ■ J FCm)  wal(n,m)  (69) 

m-O 

n ■ 0,1,2,. . . ,N-1 

where  F(m)  is  the  m th  normalized  Walsh  coefficients,  f(n)  Is  the 
discrete  Input  vector,  and  val(m,n)  Is  the  m th  Walsh  function 
(Ref  77:457).  It  should  be  noted  that  since  the  Walsh  matrix  Is 
orthogonal,  the  following  relationship  holds 

wal(n,m)  •*  wal(m,n)  (70) 

Using  matrix  notation,  the  Walsh  transform  matrix  equation  Is 
given  by 

[A]  - I [W]  [F]  I (71) 

where  [F]  is  a column  vector  of  sample  values  of  the  input  signal, 

[W]  Is  the  Walsh  matrix,  and  [A]  Is  the  column  yector  or  Walsh 
coefficients  (Ref  84:7).  \ 

Similarly,  the  "Hadamard"  transform  is  given  by 

[A]  - ^ [H]  [F]  (72) 

f^ere  [H]  Is  the  Walsh  matrix  In  natural  order  (Ref  78:178), 


43 


The  function  F(m)  or  [A],  therefore,  represents  the  "sequency" 
spectrum  of  f(n)  in  the  same  sense  that  a set  of  Fourier  coefficients 
represents  a frequency  spectrum  (Ref  44:51j52). 

The  computational  load  for  either  Equation  71  or  Equation  72  can  be 
seen  to  be  N(N-l)  additions  (multiplication  by  -1  is  not  really  a 
multiplication  but  just  a "sign"  change  and  an  addition). 

However,  Good  has  developed  a matrix  factorization  technique  which 
leads  to  a "fast"  transform  algorithm  (Ref  8:16,  17).  Good's  technique 
can  be  used  to  factor  Kroneker  matrices  such  as  in  Figure  18a  and  b. 
Matrices  of  order  N = n^  can  be  factored  into  p matrices  of  order  N 
.(Ref  8:17).  If  the  matrix  to  be  factored  has  been  generated  by  the 
Kroneker  product  of  Identical  matrices,  then  its  factors  will  also  be 

3 

identical  (Ref  84:8).  The  factors  of  a Walsh  matrix  of  order  N = 8 =2 
are  shown  in  Figure  19  (Ref  84:9).  Since  the  matrix  factors  include 
many  zero  elements,  the  number  of  computations  is  reduced  to  Nlog^N 
additions  (Ref  84:8). 

A flow  diagram  of  the  Fast  Walsh  transform  is  given  in  Figure  20 
(Ref  84:9).  Similarly,  a flow  diagram  for  the  Hadamard  transform 
is  shown  in  Figure. 21  (Ref  76:178).  The  recursive  structure  of  the 
diagrams  leads  to  an  efficient  programming  of  the  algorithm  on  a 
digital  computer  (Refs  78:179;  84:8;  1:276).  The  coefficients  of  the 

Fast  Walsh  transform  (FWT)  are  in  "bit  reflected"  order  and  those  of  the 
Fast  Hadamard  transform  (PHT'  are  in  sequency  order.  Therefore,  the  FWT 
requires  a reordering  procedure  which  adds  approximately  15  to  20%  more 
to  execution  time  (Ref  51:204).  Several  algorithms  exist  for  converting 
from  bit  reflected  order  to  sequency  order  (Ref  53:16). 
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f(0) 


Computing  Direction 


Dyadic  Convolution . The  dyadic  convolution  of  two  functions  f(t)  and 


g(t)  is  d(<flned  as 


h(t)  - f ® g 


■ J"  fCT)g(t«T)dT  (73) 

where  t denotes  dyadic  or  logical  convolution  and  t 6 t denotes  addition 
modulo-2  (Ref  41:616-617).  The  discrete  dyadic  convolution  of  two 
sequences  f(n)  and  g(n)  of  length  N Is  defined  by 

, N-1 

h(n)  f(l)g(nei),  n 

” 1-0 

(74) 

- 0,1 N-1 


Logical  Wlener-Khlntchlne  Theorem.  If  the  Walsh  transform  of 
f(n)  and  g(n)  is  defined  as  F(s)  and  G(8),  respectively,  then  the 
following  property  Is  true 


h(n)  - f 0 g 


F(8)  • G(b) 


where  s represents  sequency  (Ref  55:19).  The  tlme-sequency  domain 
logical  Wlener-Khlntchlne  theorem  Is  defined  as 


R(  ) - f e f 


F(s)-F(8) 


- F(8)- 


fdilch  Is  analogous  to  the  "arithmetic”  Wlener-Khlntchlne  theorem  of  the 
time-frequency  domain  (Refs  55:19;  72:271;  1:615). 


; 


Power  Spectral  Density.  The  sequency  power  spectrum,  P(s),  is 


easy  to  generate  from  the  Hadamard  transform  and  has  a physical 
interpretation  quite  similar  to  the  frequency  power  spectrum  (Ref  44:51). 
The  sequency  power  spectrum,  as  defined  by  Harmuth  is  given  by 


P(s) 


^2k-l  ^2k  ^ 

*Li  “ - 


,N/2-l 


(77) 


vjhere  2^  and  are  the  odd  and  even  sequency  discrete  Walsh 
transform  coefficients  (Refs  44:51;  64:93).  However,  this  spectrum  is 
not  Invariant  to  time  shifts  (Refs  64:92;  41:617).  Polge  et  al. 
state  that  the  variation  of  the  sequency  power  spectrum  with  the  time 
axis  position  of  th  input  data  is  a serious  drawback  in  signal  processing 
activities  unless  only  gross  spectral  features  are  desired  or  the 
possibility  of  time  synchronization  exists  (Ref  64:93).  Andrews 
and  Casparl,  in  their  work  on  generalized  spectral  analysis,  demonstrate 
the  "shift  variance"  of  the  Walsh  transform  relative  to  the  "shift 
invariant"  Fourier  transform  (Ref  8:24).  Figure  22a  shows  the  spectrums 
of  a block  pulse  and  Figure  22b  shows  the  spectrums  of  the  block  pulse 
shifted  relative  to  the  time  origin  (Ref  8:24). 

A time  shift  Invariant  power  spectrum  can  be  generated  by  defining 
the  power  spectrum  as  the  Hadamard  transform  of  the  autocorrelation 
function  (Equation  76)  and  noting  that  the  autocorrelation  function 
is  invariant  to  time  shifts  (Ref  64:93).  If  F(s),  s = 0,1,..., N-1, 
is  the  Hadamard  transform  of  the  sequence  f(n),  then  the  time  invariant 
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Figure  "22.  (a)  Orthogonal  Decomposition  of  a Block  Pulse 
(Fourier  to  Walsh  Transition) . (b)  Orthogonal  Decomposition^ 

of  a Shifted  Block  Pulse  (Fourier  to  Walsh  Transition) . 


power  spectrum  is 

N-1 

KAs)  - I Q.  F(8)  (79) 

1 ^ 

where  the  matrix  Q,  made  up  of  elements  Q.  , Is  dependent  on  F(s) 

' j®  . 

(Ref  62:93).  The  computation  of  the  power  spectrum  using  the 

\ 

autocorrelation  function  is  time  consuming,  and  the  high  speed 
advantage  of  the  Hadamard  transform  is  lost  (Ref  62:93). 

1 ■ ■ 

Simulation  Program  Test.  The  possibility  of  using  the  Walsh  transform 

In  the  simulation  program  TGTID  was  investigated  by  substituting  the  FFT 

subroutine  with  a fast  Walsh  transform. 


The  test  was  made  with  S broutlne  FHTl  (Appendix  D)  substituted 
for  Subroutine  FFT6  (Appendix  C) . Appropriate  :aodl float Ions  were  made 
to  Subroutine  CIMAGE  (Appendix  B,  Version  2),  and  Subroutine  ROLL 
(Appendix  B)  was  replaced  by  Subroutine  WROLL  (Appendix  B) . The 
"standard"  Walsh  power  spectrum  was  computed,  and  each  spectral 
element  of  the  resultant  Image  matrix  was  compared  with  a threshold  and 
plotted  If  It  exceeded  the  threshold.  The  threshold  was  varied  to 
determine  Its  effect  on  the  Image.  The  results  of  the  first  test  are 
shown  In  Figure  23-27,  beginning  on  page  53.  It  can  be  seen  that 
the  image  does  not  look  entirely  like  that  of  an  Image  generated  using 
a Fourier  transform.  The  effect  of  time  shifting  the  Input  was  tested 
by  "repositioning"  the  four  point  scatterers  (Appendix  H,  Data  Set  2). 

A new  set  of  images  were  constructed  and  are  shown  In  Figures  28-32, 
beginning  on  page  58 . It  can  be  seen  that  the  scatterer  "information" 
has  changed  and  that  the  two  sets  of  images  have  only  little  similarity. 
The  second  data  set  was  used  with  the  FFT  subroutine  reinserted  Into 
the  simulation  program.  It  ca  be  seen  that  there  Is  no  change  In  the 
image,  except  in  its  location  in  the  "viewing  field",  and  the  result 
Is  shown  In  Figure  33,  on  page  63  • 

This  difficulty,  as  Blachman  observes,  is  attributed  to  the  fact 
that  after  time  shifting,  a Walsh  function  generally  becomes  the  sum 
of  an  Infinite  number  of  Walsh  functions  while  a sinusoid  simply 
turns  Into  the  sum  of  a sine  and  cosine  of  the  same  frequency.  Thus 
a change  In  time  scale,  or  a shift  of  the  time  origin  usually  will 
grossly  alter  a Walsh  spectrum  but  has  no  effect  on  the  Fourier  spectrum 
(Ref  16:347). 
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It  Is  concluded  that  since  the  Valsh  power  spectrum  computed  by  the 
direct  method  Is  not  time  shift  Invariant  and  since  no  fast  algorithm 
exists  for  the  second  method  and  since  there  Is  little  possibility 
of  time  synchronization,  the  Walsh  transform  can  not  be  utilized 
in  this  application. 
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IV.  Walsh  Domain  to  Fourier  Domain  Conversion 


A technique  developed  by  Andrews  and  Casparl  (Ref  9)  implements 
a fast  Fourier  transform,  a fast  Hadamard  transform,  and  a variety 
of  other  orthogonal  decompositions  suggesting  a generalized  spectral 
analysis.  Their  results  imply  that  a relationship  exists  between 
the  Walsh  domain  and  the  Fourier  domain.  Robinson  (Refs  72,  73) 
presented  a derivation  of  a recursive  relationship  between  the  arith- 
metic and  logical  autocorrelation  functions  of  a wide  sense  stationary 
process.  This  work  was  based  on  a theorem  discovered  by  Gibbs  and 
proved  by  Plchler  (Ref  38).  Siemens  and  Kltal  (Refs  79,  80)  and 
Blachman  (Refs  15,  16)  describe  schemes  for  converting  from  Walsh 
coefficients  to  Fourier  coefficients.  Both  approaches  seem  promising 
because  of  the  computational  speed  advantage  of  the  FWT/FHT  as  compared 
to  the  FFT  and  will  be  analyzed  In  this  section.  However,  It  is  shown 
that  Robinson's  approach  Is  Incorrect  and  the  second  approach  not 
feasible  for  use  In  radar  target  Identification. 

Walsh  Power  Spectrum  to  Fourier  Power  Spectrum 

Robinson  defines  the  Walsh  power  spectrum  of  a sequence  of  rt.nd(^ 
samples  as  the  Walsh  transform  of  the  "logical"  autocorrelation 
function  of  the  random  sequence,  where  the  logical  autocorrelation 
function  Is  defined  In  a similar  form  as  the  "arithmetic" 
tlon  function  (Ref  72:271).  Robinson  asserts  that  the  Fourier  power 
spectrum,  which  Is  defined  as  the  Fourier  transform  of  the  arithmetic 
autocorrelation  function,  can  be  obtained  from  the  Walsh  power  spectrum 
by  a linear  transformation  (Ref  72:271).  The  chain  of  transformations 
can  be  summarized  as 
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I 


Fourier 


Arithmetic 


Tower 


Spectrum 


Walsh 


Power 


I Spectrum 


(Autocorrelation) 


Function 


Ta-l  (79) 


Logical 

Autocorrelation 


Ftinctlon 


Hcvever,  the  procedure  to  find  the  transformation  matrix,  T,  was 
found  to  be  Incorrect.  Robinson  defines  the  logical  autocorrelation 
function  as 


L(“>(k)  - I i x(j  •k)x(j) 
J-0 


k - 0,1,2,. 


where  x(J),  J ••  0,1,2, ... ,N-1,  Is  a random  sequence  of  length  N » 2° 
and  represents  a window  or  block  of  N samples  of  discrete  random  pro- 
cess. The  logical  autocorrelation  function  Is  then  defined  as  the 
expected  value  of  the  local  logical  autocorrelation  function  of 
Equation  80 

L(k)  E{L(“>(k)}  (81) 


where  the  expectation  operator  E denotes  the  ensemble  average  of  M 
local  logical  autocorrelation  functions  (Ref  73:299) 


m-l 

k - 0,1,2, ...,N-1 
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Robinson  then  defines  the  arithmetic  autocorrelation  function  as 
as  even  function  of  time  difference  only 


E{x(j+k*)x(j)}  - R(k*) 


(83) 


where  k*  Is  the  time  shift  (Ref  72:272). 

If  Equation  81  Is  written  In  matrix  form  using  the  Indices  of 
Table  I,  Robinson  asserts  that  a linear  combination  of  R(k)  Is  obtained 
for  each  L(k'.  As  an  example,  Robinson  presents,  for  N ■ 4,  the 
following 


1(0) 

x(0) 

x(l) 

x(2) 

x(3) 

x(0) 

L(l) 

1 

x(l) 

X(0) 

x(3) 

x(2) 

x(l) 

- E < 

A 

L(2) 

4 

x(2) 

x(3) 

x(0) 

x(l) 

x(2) 

un 

w 

x(3) 

x(2) 

x(l) 

x(0) 

x(3) 

(84) 


Robinson  shows  that  the  first  row  which  corresponds  to  L(0)  yields  the 
correlation  of  N samples  for  zero  time  shift,  thus  L(0)  R(0),  which 
Is  true.  Robinson  also  states  that  L(l)  ■>  R(l)  which  Is  not  true. 

The  standard  definition  of  the  arithmetic  autocorrelation  function  Is 
given  as  (Ref  88:13) 


R<t)  - ^ I x(n)x(n+|T|) 
" n-0 


(85) 


Using  Robinson's  relation.  Equation  84,  L(l)  Is  found  to  be 

L(l)  - l/2[x(l)x(0)  + x(3)x(2)]  (86) 

Using  Equation  85,  It  can  be  seen  that 
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(87) 


R(l)  - l/4[x(0)x(l)  - x(l)x(2)  + x(2)x(3)  + x(3)x(0)] 


Table  I 


Bitwise  Modulo  2 Addition,  Given  by  j ^ k,  for 
Integers  Between  0 and  3 


j 

k 

0 

1 


0 12  3 

0 12  3 

10  2 3 


2 2 3 0 1 

3 3 2 1 0 

(Ref  72:272) 


A J 
'-i 


Table  II 


Time  Shift  k*.  Given  by  k*  - (J  • k)  - j,  for 

Integers  j and  k Between  0 and  3 

i ■' 

» i 

■j 

J 

0 

1 

2 3 

k 

0 

0 

0 

0 0 

1 

1 

-1 

1 -1 

■’  / 

2 

2 

2 

-2  -2 

/■ 

3 

3 

-1 

-3  -3 

/ .. 

(Ref  72:272^ 

therefore. 


L(l)  ^ R(l) 


(88) 


Robinson's  error  arises  from  the  confusion  of  the  use  of  k*  in 


\ 

\ 
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Equation  83  as  a dummy  variable  and  the  use  of  k*  In  Table  II  as  the 
time  shift  for  computing  R(k).  This  Is  an  Incorrect  use  of  the  defi- 
nition of  the  arithmetic  autocorrelation  function  given  In  Equation  85, 
that  Is,  the  time  shift  cannot  be  varied  when  computing  the  arithmetic 
autocorrelation  function  for  a particular  time  shift.  Robinson  uses 
Table  II  to  find  a specific  k*  for  each  J In  the  summation  In  com- 
puting R(k}. 

It  Is  concluded  that  this  approach  is  Incorrect  and  that  there 
Is  no  linear  relationship  between  the  Walsh  power  spectrum  and  the 
Fourier  power  spectrum  of  a signal. 

Walsh  Series  to  Fourier  Series  Coefflcelnts 

Siemens  and  Kltal,  and  Blachman  have  shown  that  the  coefficients 
of  the  Walsh  series  of  a function  can  be  used  to  derive  the  correspon- 
ding Fourier  series  coefficients.  The  conversion  equation  for  each 
Fourier  coefficient  is  In  the  form  of  an  Infinite  summation  of 
products  of  constants  and  the  Walsh  coefficients  (Ref  79;295).  They 
assume  that  the  signal  Is  frequency-limited  so  that  precise  evaluation 
of  the  Fourier  coefficients  In  terms  of  Walsh  coefficients  is  possible 
and  that  the  highest  normalized  frequency  component  (harmonic)  N and 
the  highest  normalized  sequency  component  M are  equal.  Thus  a finite 
number  of  Walsh  coefficients  can  be  used  for  computation  of  the  Fourier 
coefficients.  The  conversion  computation  Is  further  reduced  if  H is 
a power  of  two  (Ref  79:295). 

Let  a function  f(0)  be  represented  by  a sequency-ordered  Walsh 
series 
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(89) 


f(e)  - Ap  + I [Aj,  cal(m,e)  + sal(m,0)] 

1B*0 

The  coefficients  Aq  end  of  the  even  and  odd  Walsh  functloris,  res- 
pectively, are  defined  by 


K - 

f f(9)cal(m,9)d6 
* 

(90) 

u 

A 

Bm- 

f(9)sal(m,9)d9 

(91) 

0 


The  same  function  f(n)  has  the  corresponding  Fourier  series 


f (fl)  ■ 2 ^ tanCOs2im9  + bQsln2Trn6l 


(92) 


n«l 


where 


a-  ■ 2 I f(6)  cos  2im6  d6 

h 

.1 


-2  I 

•'0 


f(9)  slu  2im6  d9 


(93) 


(94) 


The  objective  Is  to  use  the  Walsh  coefficients  A^  end  B^^  to  derive  a^ 

and  b_ . 1 

“ ! 

Siemens  and  Kltal  claim  that  the  even  terms,  a„,  of  the  Fourier 

n 

series  of  a signal  are  functions  only  of^^the  even  terms,  Ajj,,  of  the 
corresponding  Walsh  series  (Ref  79:295).  1 Similarly,  b^  terms  depend 
only  on  B^,  terms.  The  even,  real  terms  ajre  derived  below;  similar 
derivations  apply  for  the  odd  terms. 

The  Walsh  to  Fourier  series  conversion  Is  derived  by  equating 
the  terns  of  each  series  (Ref  79:296) 
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I aj^cos2Trn0  = I A^jCal(m,0) 
n=l  m=l 


(95) 


The  cal  functions  are  expanded  into  sets  of  equivalent  Fourier  series 
expressions  whose  terms  have  coefficients  a^  ^ where 


a =2 
n,m 


cal(m,0)cos2'iTn0d0 


(96) 


If  a represents  the  nx]  matrix  of  the  set  {a^}  as  n ->■  and  A repre- 
sents the  mxl  matrix  of  the  set  {A^^}  as  m then 


a = F A 


(90 


If  only  a finite  number  M of  Walsh  coefficients  are  known,  then  a^  can 

A 

only  be  approximated  as  a^^,  where 


M 


m=l 


(99) 


is  the  n th  Fourier  coefficient  of  cal  (ra,0).  The  mxn  matrix  of  the 
set  {a^  jjj)  is  denoted  F^  (Ref  79:296).  In  the  expansion  of  the  right 
hand  side  of  Equation  96,  terms  containing  cos2iin0  are  grouped, 
yeilding  a^  values  given  by 


00 


(97) 


However,  if  the  function  f(n)  is  either  frequency-limited  or  sequency- 
limited , than 


^n  ^ 


(100) 
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m: 


The  coefficients  a can  be  considered  the  Fourier  coefficients  of  tb& 
n 

frequency-limited  or  sequency-llmlted  function  (Ref  79:296). 

T 

The  next  step  is  to  find  the  set  of  F , the  Fourier  co- 

efficient' set  of  cal(m,6}.  Siemens  and  Kltal  developed  an  alternative 
expression  for  the  Fourier  transform  of  a Walsh  function  which  differed 
from  earlier  expressions  In  that  It  Incorporated  the  Gray  code  repre- 
sentation of  the  order  of  the  function  (Refs  79:81;  15:349). 
Additionally,  the  expression  le  nonrecurslve.  llie  definition  of  the 
Fourier  transform  of  a Walsh  function,  val(m,6).  Is  defined  as 


r[wal(m,6)]  ■ [ wal(m,0)exp(j2rfe)d6 

Jo 


(101) 


The  even  and  odd  Walsh  functions,  cal(m,6)  and  8al(u,6),  respectively, 
have  the  Fourier  transforms 


Flcal(8,0)]  - C(f,8) 


■ f cal(8,6)co82rf0d0 

» ft 


(102) 


F[8al(8,0)l  - JS(f,8) 


A 

> I 8al(8,0)8ln2irf0d0 

Jo 


(103) 


where  f and  s are  normalized  frequency  and  sequency.  Since  the  Waleh 
functions  are  discontinuous,  evaluation  of  Equation  101  - Equation  103 
would  normally  involve  a summation  of  integrals  (Ref  80:81). 

Siemens  and  Kltal  state  that  It  is  convenient  to  view  a contln- 


-t  :v 


.-■iV' 


M - 


uous  Walsh  function  as  the  convolution  of  the  sequence  of  unit  Ish* 
pulses  that  form  the  discrete  Walsh  function  over  the  Interval 
- 0 < 6 1 with  a rectangular  pulse  of  unit  magnitude  and  the  width 

of  1/2"  equal  to  the  spacing  of  the  unit  Impulse,  where  M 1^  the  number 
of  bits  In  the  binary  representation  of  m.  The  Fourier  transform  of 
the  Walsh  function  Is  then  the  product  of  the  trensforms  of  the  discrete 
Walsh  function  and  the  rectangular  pulae  (Ref  80:81).  This  relation 
Is  then 

Flwal(m,0)]  - (-l)®0(_i)0.  [“n^08(5i  - g^  |)  ] •8lnc(f/2**)  (104) 


where  g^  Is  the  bit  In  the  Gray  code  representation  of  m,  and  a Is 
the  number  of  Gray  code  bits  of  value  ONE,  and 


sinc(f/2M)  A 8ln(  f/2”)/(  f/2**) 


(105) 


If  the  cal  and  sal  functions  are  given,  the  order  m Is  found  from 


wal(m,6)  * cal(s,6), 
in  " 2s 


(106) 


wal(m,0)  ■ 8al(8,6), 
m • 2s  - 1 


(107) 


Equation  99  Is  sufficient  to  compute  the  Fourier  coefficients  of 
f(n)  assuming  the  set  has  b^en  found  and  Equation  104  has  been 

used  to  determine  the  F^’s  necessary  to  find  Fourier  coefficients.  A 
total  of  2M  Walsh  coefficients  are  necessary  to  find  the  complex 
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Fourier  coefficients.  Therefore,  to  find  A requires  2Mlog2M  real 
additions  and  subtractions.  A total  of  2M^  storage  locations  are 
required  for  the  Fourier  coefficient  sets  of  cal(m,6)  and  sal(m,e). 

The  computational  load  to  find  the  Fourier  coefficients  of  f(n) 
can  be  seen  to  be  2M^  multiplies  and  (2M^  + 2mLog2M)  additions.  This 
Is  compared  to  4MLog2M  real  multiplies  plus  2I{H  additions  to  compute 
the  fast  Fourier  transform. 

It  can  be  seen  that  the  Semes  conversion  approach  requires 
greater  computational  effort  than  the  direct  fast  Fourier  transform 
method  and  a very  large  storage  requirement  for  the  transformation 
matrices  (F^).  It  Is  therefore  concluded  that  this  approach  is  not 
feasible  for  Implementation  In  this  problem. 
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V.  Comparison  of  FFT  Algorithms 

It  is  seen  in  Section  IV  that  many  orthogonal  transforms  are  not 
feasible  for  use  in  tactical  target  identification  signal  processing. 
This  section  provides  a brief  description  of  the  fast  Fourier  transform 
(FFT)  and  then  selects  the  optimal  implementation  of  it  for  this  appli- 
cation. Performance  criteria,  based  on  execution  time,  memory 
requirements,  and  error  size,  are  developed  and  serve  as  a basis  of 
comparison  among  the  implementations  considered.  Two  general  t}rpes  of 
implementations,  floating  point  and  fixed-point,  are  presented,  tested, 
and  evaluated  using  the  criteria. 

Several  different  "floating  point"  FFT  subroutines  are  tested  using 
the  simulation  program  TGTID  as  the  "driver"  program  to  determine  their 
relative  utility.  Two  of  these  FFT  subroutines  are  modified  to 
execute  more  efficiently  then  one  FFT  is  selected  based  on  the  criteria 
as  the  most  feasible  in  the  simulation  program.  The  FORTRAN  IV  code  of 
the  FFT  subroutines  is  listed  in  Appendix  C. 

A "fixed-point"  FFT  implementation  is  derived  and  demonstrated  to 
be  feasible  in  this  application.  It  is  considered  to  be  the  optimal 
approach  for  this  problem.  Several  variations  of  the  fixed  point 
FFT  are  Implemented  with  the  FORTRAN  IV  listings  given  in  Appendix  E. 
One  fixed-point  FFT  subroutine  is  selected  as  the  best  choice  for  use 
in  the  simulation  program. 

Background 

The  FFT  algorithm  is  a highly  efficient  procedure  for  computing 
the  Discrete  Fourier  transform  CDFT) 


(108) 


II-l 

X(k)  = I x(n)exp(-j27tnk/N) 
n=0 

for  k •=  0,1 ,2, . • . ,N-1,  where  x(n),  n 0,1,2,. ..  ,N-1  Is  an  equally 
spaced  complex  data  series  (Ref  62:100).  It  takes  advantaj^e  of  the  fact 
that  the  calculation  of  the  Fourier  coefficients  can  be  carried  out 
Iteratively,  which  results  in  a considerable  savings  of  computational 
time  (Ref  23:45).  The  DFT  is  a reversible  mapping  or  unitary  operation 
for  time  series  and  has  the  same  mathematical  properties  as  those  of 
the  Fourier  integral  transform  given  in  Equation  8 (Ref  23:45).  In 
particular,  it  defines  a frequency  spectrum  of  a time  series. 

If  digital  analysis  techniques  are  used  for  analyzing  a continuous 
waveform  then  it  is  necessary  that  the  data  be  sampled  at  equal  time 
intervals  in  order  to  produce  a time  series  of  discrete  samples  which  can 
be  used  in  a digital  computer.  So  that  the  time  series  completely 
represents  the  continuous  waveform,  it  is  necessary  that  the  waveform  be 
band-limited  and  sampled  at  twice  the  highest  frequency  in  the  waveform 
(Ref  24:45).  This  rate  is  known  as  the  Nyqulst  rate  and  the  equispaced 
samples  are  known  as  Nyqulst  samples. 

The  speed  of  the  FFT  algorithm  depends  on  the  factorabllity  of  N 

M 

N - n n (109) 

1 » 1 ^ 

and  then  decomposing  the  transform  into  M stages  with  N/n^  transformations 
of  size  n^  within  each  stage  (Ref  18:93).  The  FFT  algorithm  is  most 
efficient  when  N = 2^  or  n = 4^/2  since  a cotnpu.':atlonal  advantage  is 
gained  by  decomposing  the  computation  of  the  DFT  into  successively 
smaller  DFT's  (Refs  54:2;  62:285-287).  The  improvement  in  computational 

efficiency  of  the  FFT  is  proportional  to  N log^N  as  opposed  to  the 
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computational  load  of  N for  the  straight  DFT  approach  (Ref  60:287). 

There  are  two  basic  FFT  algorithms.  The  first,  called  "Decimntlon- 
in-tlme",  derives  its  name  from  the  fact  that  in  the  process  of 
arranging  the  computation  into  smaller  transformations,  the  sequence 
x(n)  is  decomposed  into  successively  smaller  sequences.  In  the  second 
general  class  of  algorithm,  the  sequence  of  Fourier  transform  coefficients 
of  X(k)  is  decomposed  into  smaller  subsequences,  hence  the  name 
"declmatlon-in-frequency"  (Ref  62:286-287).  When  either  one  of  the 
two  algorithms  is  used,  the  ordering  of  the  sequence  (either  input  or 
output)  must  be  taken  into  account.  If  the  declmlnatlon-ln-frequency 
algorithm  is  used  the  output  coefficients  will  be  in  bit-reversed 
order;  and  if  the  declmatlon-ln-tlme  algorithm  is  used,  the  input  samples 
must  be  put  into  bit-reversed  order.  Therefore,  in  addition  to 
performing  the  FFT,  a reordering  procedure  must  also  be  implemented  which 
accounts  for  approximately  15-20Z  of  the  execution  time  of  the  FFT 
process  (Ref  64:17).  Figure  34  and  Figure  35  show  the  signal  flow 
graphes  of  the  declmatlon-in-tlme  and  decimatlon-ln-frequency  algorithms, 
respectively. 

^ The  literature  is  rich  in  material  on  the  FFT:  its  derivation, 

\ 

ltd  properties,  its  uses,  and  its  pitfalls.  The  reader  is  referred 

particularly  to  Cooley  and  Tukey  (Ref  25).  Bergland  (Ref  12), 

\ 

I 

Cochran  et  al  (Ref  24),  Gentleman  and  Sande  (Ref  36),  Singleton  (Refs  81, 

82), I and  Oppenheim  and  Schafer  (Ref  62).  The  Bibliography  of  this 

thesis  contains  several  other  related  and  excellent  references. 

Criteria 

0 

The  criteria  generally  applied  to  the  evaluation  of  an  FFT  algorithm 
Include  program  execution  time,  accuracy,  storage  requirement  , and 
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■Igure  34.  Flow  Graph  of  Complete  Decimatlon-ln-tlme  Decomposition  of  an  Eight-point  DFT  Computation 


adaptlblllty  to  the  digital  computer  to  be  used  (Ref  34:1)  Bergland 
points  out  that  the  number  of  variations  of  the  FFT  algorithm  appears 
to  be  directly  proportional  to  the  number  of  people  using  It  and  that 
most  of  these  Implementations  are  based  on  either  the  Cooley-Tukey  or 
the  Sdnde-Tukcy  methods,  but  are  formulated  to  exploit  properties  of 
the  time  series  being  analyzed  or  the  properties  of  the  computer  being 
used  (Ref  12:50).  Ferrle  states  that  best  accuracy  Is  achieved  only 
at  the  expense  of  Increased  execution  time  and  storage  (Ref  34:10). 

In  essence,  then,  no  single  FFT  algorithm  represents  a "best"  choice:  | 
it  depends  upon  the  application  (Ref  71:25). 

For  this  application,  short  execution  time  Is  considered  the  most 
liq>ortaitt  criterion,  and  the  next  most  Important  Is  small  storage 
requirements.  Since,  only  major  spectral  lines  are  of  Interest,  then 
only  coarse  calculations  are  Important,  therefore,  accuracy  is  not  | 

considered  an  Important  criterion.  Additionally,  the  type  of  digital 
computer  for  Implementation  is  the  DEC  FDF-11/40,  and  Its  high  speed 
Integer  arithmetic  capability  motivated  the  development  of  a fixed- 
point  FFT  subroutine. 

Floating  Folnt  FFT  Algorithm  Comparison 

The  eight  FFT  subroutines  evaliiated  are  briefly  described  In  this 
section.  A.1^  the  subprograms  are  written  In  FORTRAN  IV.  In  their 
original  form  most  were  named  "FFT";  for  this  thesis  they  are  designated 
FFTl  through  FFT6.  Two  of  the  subroutines,  FFTl  and  FFT4,  are  modified 
to  improve  their  efficiency.  All  the  subroutines,  except  FFT2,  perform 
the  transform  In-place.  And  all  handle  complex  exponentials  In  terms  of 
sine  and  cosine  functions  according  to  Euler's  rule 

exp(-j0)  ■ cos e -jsln  6 (110) 
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FFTIA  and  B (Ref  62:331).  FFTlA  and  B are  radix-2,  decimation- 
In-tlme  transform  subprograms.  They  can  handle  N equal  a power  of  two 
points.  The  data  are  treated  as  complex,  and  the  trigonometric 
functions  are  computed  as  they  are  used.  They  are  simply  implemented 
and  straightforward.  Computations  are  performed  in-place  and  calcu- 
lations use  complex  arithmetic.  Memory  storage  requlremento  depend 
only  on  N.  Initial  bit  reverse  ordering  is  done  within  the  subroutine. 

FFTIB  is  made  more  efficient  by  performing  the  first  state 
separately  from  the  rest  of  the  stages  by  exploiting  the  fact  that  the 
vrlue  of  ••  exp(-j2irnk/N)  is  zero  for  the  first  stage  butterflies, 

the  basic  unit  of  computation.  Therefore,  the  cosine  and  sine  values 
are  one  and  zero,  respectively,  and  the  arithmetic  operations  are 
simple  additions  and  subtractions.  Figure  34  shows  this  property. 

FFT2  (Ref  71:21).  FFT2  is  a radix-2  decimation-in-tlme  algorithm. 
It  uses  two  arrays  to  perform  the  transform,  so  that  no  separate 
reordering  routine  is  required,  but  the  storage  requirement  is 
doubled.  Complex  arithmetic  is  used.  It  is  known  for  its  accuracy 

but  its  execution  is  slow  (Ref  71:25). 

FFT3  (Ref  71:21).  FFT3  is  a simple,  straight  forward  radlx-2 
implementation  of  the  decimatlon-ln-frequency  algorithm.  It  uses 
real  arithmetic  to  compute  the  butterfly.  It  requires  Subroutine 
RBITS  to  reorder  the  bit-reversed  Fourier  coefficients.  Computations 
are  performed  in-place.  Internal  storage  requirements  for  this 
subroutine  are  greater  than  those  of  the  previous  subroutines. 

FFT4A  and  B (Ref  62:332).  FFT4A  performs  the  transform  approx- 
imately the  same  as  FFT3  except  that  the  indexing  scheme  is  simpler. 
Complex  arithmetic  is  used  and  calculations  are  performed  in-place. 
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Reordering  is  accomplished  within  the  subroutine.  FFT4B  Is  modified 
to  make  It  more  efficient  by  performing  the  last  stage  separately, 
similar  In  concept  to  FFTIB. 

FFT5  ' (Ref  52:27).  FFT5  Is  based  on  the  original  Cooley-Tukey 
FORTRAN  subroutine.  It  Is  a radlx-2  transform  subprogram  and  performs 
computations  In-place.  It  Is  neither  a fast  or  efficient  Implementa- 
tion of  the  FFT  algorithm  (Ref  54:28). 

FFT6  (Ref  45:A-3).  FFT6  Is  a mixed-radix  declmatlon-ln-frequency 
algorithm.  It  can  handle  2^^  points  which  can  be  Increased  with  minor 
modification  to  the  subroutine.  It  requires  significantly  more  Inter- 
nal storage  than  any  of  the  other  FFT  subroutines,  but  It  Is  also 
known  for  Its  speed  (Ref  82). 

Evaluation  of  the  Floating  Point  FFT  Subroutines 

The  eight  FFT  subroutines  were  tested  In  the  simulation  program 
TGTID  with  appropriate  modifications  made  to  Subroutine  CIMAGE 
(Appendix  B). 

Execution  times,  memory  requirements,  and  the  number  of  FORTRAN  IV 
statements  for  each  FFT  subroutine  are  summarized  In  Table  III.  The 
number  of  FORTRAN  IV  statements  Is  the  total  number  of  statements 
required  to  perform  the  transform  and  reordering  routines,  but  does 
not  Include  comment  cards.  Storage  requirements  were  obtained  from 
the  cross-reference  maps  produced  by  the  CYBER  SCOPE  Operating  System. 
Execution  time  was  determined  by  using  the  CEC  library  subroutine 
SECOND  (CP)  which  returns  time  In  seconds  to  three  decimal  places.  In 
Subroutine  CIMAGE,  the  FFT  subroutine  Is  called  256  times  (NTM  ■ 2 
and  MBSZ  * 128;  It  Is  called  128  times  for  NTM  » 1 and  128  times  for 
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NTM  = 2).  The  FFT  eize  is  128  points  (IJBSZ  = 128  = 2^20^°,  where 
NBORD  = 7) . Thus  each  FFT  tested  performed  256  128-point  transforms 
to  produce  the  image.  The  execution  tine  was  measured  by  taking  a 
time  "hack"  (STIME)  just  before  the  FFT  subroutine  was  called  and  just 
after  it  returned  (ETDIE)  to  the  calling  subprogram.  The  difference 
between  the  two  tine  hacks  was  computed  (RTIME  = ETIME  - STIME)  and 
added  to  the  total  accumulated  execution  time  (TTIME) . The  final 
value  of  TTIME  was  the  total  execution  time  to  compute  the  256  128- 
point  FFT  subroutines  needed  to  generate  one  image.  The  execution 
time  entry  in  Table  III  is  an  average  of  the  total  execution  time  for 
each  FFT  subroutine  to  generate  at  least  three  separate  images. 

Results  of  Floating  Point  FFT  Comparison 

From  Table  III,  it  can  be  seen  that  FFTIB  and  FFT  6 are  the 
fastest  subroutines.  FFTIA  required  the  least  amount  of  memory. 
However,  it  is  concluded  that  FFTlB  is  the  best  floating  point  FFT 
subroutine  based  on  its  speed  and  storage  requirements.  Figures 
36-39  show  the  output  inages  generated  by  FFTlB  for  increasing 
threshold  settings.  Images  generated  using  FFT6  are  shown  in 
Figures  9 through  lA. 

Fixed-Point  FFT 

The  implementation  and  use  of  a fixed-point  FFT  is  motivated  by 
several  factors.  The  DEC  PDP-11 /AO  minicomputer  to  be  used  in  the 
actual  applications  work  of  radar  target  imaging  processes  an 
extremely  fast  integer  arithmetic  proce.ssor  and  has  no  floating 
point  processor.  Tlie  16-blt  word  is  considered  sufficient  to  provide 
an  adequate  resolution  since  only  gross,  relative  spectral  line 
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Figure  38, 


Subroutine  FFTIB;  Threshold  = 1100 


magnitudes  are  really  required  for  this  application.  Fixed  point  division 
by  two  is  easily  performed  by  shifting  the  binary  point  to  the  left. 

The  fixed-point  FFT  is  based  on  the  16-bit  word  of  the  DEC  PD?-il/40. 
The  data  word  is  divided  into  two  parts:  the  lower  15  bits  for  the 
numbers  are  scaled  so  that  the  binary  point  lies  at  the  extreme  left. 

The  input  data  must  be  put  into  the  proper  format.  First,  each 
point  is  multiplied  by  2^^  - 1 ««  32767  to  put  it  into  a scaled  integer 
format.  Second,  each  point  is  divided  by  the  total  number  of  points, 

N,  in  the  transform  (in  this  case  N *■  128)  to  negate  the  gain  that 
is  generated  as  the  computation  of  the  FFT  progresses  from  stage  to 
stage  to  prevent  the  possibility  of  an  overflow  when  two  points  are 
added  together  in  the  butterfly  calculation.  The  real  and  Imaginary 
parts  of  each  data  set  are  scaled  and  are  sorted  in  two  separate,  real 
and  Imaginary,  arrays,  respectively. 

The  trigonometric  functions  are  computed  and  then  multiplied  by 

32767  to  integer  scale  each  sine  and  cosine  value.  When  the 

trigonometric  terms  are  multiplied  by  the  difference  of  two  points  in 

the  computation  of  a butterfly,  the  product  results  in  a 31-blt 

number  plus  sign,  which  is  greater  than  one.  To  rescale  this  number, 

15 

it  is  divided  by  2 = 32768  or  more  simply,  the  binary  point  is 

shifted  15  jisces  to  the  left  and  the  16  least  significant  bits 
are  truncated.  This,  likewise,  prevents  an  overflow  condition  from 
occurring.  The  basic  butterfly  computational  algorithm  (decimation- in- 
frequency) is  shown  in  Figure  40,  and  the  corresponding  equations  are 


T1  - X (p) 
a ^ 

- 

(111) 

T2  - Y^(p) 

(112) 
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Figure  40.  Flow  Graph  of  Two-Point  Integer  FFT 
Butterfly  Using  Real  Arithmetic 


Vl*''  ■ (113) 

VlW  ■ * ^(3)  (IM) 

Xji^lCq)  * (cos(2  r)*32767*Tl  + sin(2  r) *32767) /32768  (115) 

Yj^l(q)  - (cos (2  r)*32767*T2  + sin (2  r)*32757)/32768  (116) 


where  X and  Y represent  the  real  and  imaginary  parts,  respectively. 

The  accuracy  of  the  fixed-point  power  of  two  FFT  algorithm  is 
addressed  by  Welch  (Ref  91) . His  error  analysis  led  to  approximate 
upper  and  lower  bounds  on  the  root-mean-square  error.  Based  on 
Welch's  findings,  it  was  determined  that  the  theoretical  upper  bound 
for  this  application  (B  = 16  and  N = 128)  is  (Ref  91:156) 
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RMS (error) 
RMS  (result) 


(117) 


*»  l.SxlQ-^ 

This  bound  Is  considered  to  be  more  than  adequate  for  radar  target 
imaging.  Welch  further  states  that  if  the  transform  is  taken  to 
•estimate  the  power  spectrum  of  a signal,  averaging  over  frequency  in  a 
single  perlodogram,  or  averaging  over  time  in  a sequence  of  perlodograms, 
or  using  simple  wleghting  will  decrease  the  error  (Ref  91:157). 

It  is  concluded  that  the  use  of  the  fixed-point  FFT  algorithm  is 
feasible  for  implementation  in  radar  target  imaging  signal  processing. 

Five  "simulated”  fixed-point  FFT  subroutines  were  written  in 
FORTRAN  IV  for  this  thesis.  They  are  based  on  the  floating  point  FFI 
algorithms  previously  discussed,  exepet  Subroutine  FFTIA.  They  are 
written  in  stimulated  form  for  use  on  the  CDC  6600  CYBER  74  System 
Instead  of  the  minicomputer  for  which  they  are  Intended.  Shifting  of 
the  binary  point  is  accomplished  with  the  use  of  the  CDC  Intrinsic 

i 

library  function  SHIFT  (A,-N) . ' 

FFTIi.  FFTIl  is  based  on  FFT.  | It  is  converted  from  complex, 
floating  point  arithmetic  to  real,  fixed-point  arithmetic.  Indexing 
and  reordering  are  not  changed. 

FFTI2.  FFTI2  is  based  on  FFTl.  It  is  not  as  efficient  to 
Implement  as  the  decimatlon-ln-frequency  algorithm. 

FFTI3.  FFTI3  is  based  on  FFT4.  Reordering  is  performed  within 
the  subroutine. 

FFTIA . (Ref  34:15-18)  FFTI4  is  based  in  part  on  an  algorithm 
presented  by  Fisher  (Ref  35) . Trigonometric  functions  are  precomputed, 
scaled,  and  stored  in  a table.  The  subroutine  is  called  once  by  the 
calling  program  and  passes  the  table  to  the  main  FFTI4  subroutine. 


\ 
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Subroutine  COSINE  (Appendix  10  is  used  to  generate  the  trigonometric 
table.  It  is  called  by  Subroutine  CIMAGE  (Appendix  E,  Version  3), 
and  used  in  Subroutine  FFT14.  An  address  loopup  routine  is  used 
to  obtain  the  proper  sine/cosine  value  from  the  table.  Subroutine 
UNSCRl  (Appendix  J)  is  called  by  FFTI4  to  reorder  the  Fourier 
coefficients.  A complete  derivation  and  floating  point  implementation 
is  reported  by  Fisher  (Ref  34). 

FFTI5.  FFTI5  incorporates  the  best  features  of  FFT4B  and  FFTI4. 

A cosine  table  is  generated  from  zero  t.^  two  pi.  FFT15  uses  a 
lookup  address  routine  to  obtain  the  required  trigonometric  values 
for  the  butterfly  computation.  The  last  stage  of  the  transform  is 
performed  separately  to  eliminate  the  lookup  and  unnecessary  multiplies. 
Subroutine  UNSCRl  is  called  to  reorder  the  fixed-point  Fourier 
coefficients. 

Evaluation  of  the  Fixed-point  FFT  Subroutines 

The  five  "simulated"  fixed-point  FFT  subroutines  were  tested  in  the 
simulation  program  TGTID  with  appropriate  modifications  made  to 
Sburoutlne  CIHAGE  (Appendix  B,  Version  3).  They  were  evaluated  in 
the  same  manner  as  the  floating  point  FFT  subroutines.  Table  IV 
summarizes  the  ntimber  of  FORTRAN  IV  statements,  the  amount  of  memory 
required,  and  the  average  execution  time  (over  at  least  three  runs  of 
program  TGTID)  for  each  fixed  point  FFT  subroutines. 

The  execution  time  was  measured  in  exactly  the  same  manner  as  the 
floating  point  FIT  subroutines.  The  first  time  "hack"  CSTIME)  was  taken 
just  after  the  data  were  converted  to  fixed-point  format  and  just 
after  the  data  were  converted  to  fixed-point  format  and  just  before  the 
fixed-point  FFT  subroutine  was  called.  The  second  time  hack  (ETIME) 
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was  taken  just  after  the  return  from  the  subroutine.  The  difference 
was  computed  (RTIME)  and  added  to  the  accumulated  total  execution  time 
TTIME)  of  the  fixed-point  FIT  subroutine.  This  time  Included  obth 
the  transform  and  reordering  routines.  It  does  not  Include  the  time  to 
conqjute  the  trigonometric  functions  of  Subroutines  FFTI4  and  FFTI5,  which 
are  only  done  once. 

Results  of  the  Fixed-point  FFT  Subroutine  Comparison 

It  can  be  seen  in  Table  IV,  that  FFTI5  is  the  fastest  fixed-point 
FFT  subroutine.  FFTI3  requires  the  least  amount  of  memory,  while  FTTH 
reuqlres  the  most  amount  of  storage.  Based  on  the  speed  and  moderate 
storage  requirements,  FFTI5  is  considered  to  be  the  best  fixed-point  FFT 
subroutine  to  use  in  tactical  radar  target  Imaging  using  the  DEC  PDP-11/40 
minicomputer.  Figures  41  through  44  show  the  image  created  using  FFTI5  for 
Increasing  threshold.  It  should  be  noted  that  the  threshold  is  set  too 
low  in  Figure  41  and  too  high  in  Figure  44. 


Table  IV 


S urinary  of 

Fixed  Point  FFT 

Performance  Statist! 

cs 

FFT 

Nuisber  of 
FORTRAN 
Statements 

Number  of 

T ^ ^ • II  j Aux  Array 

Instruction  Words 

(Total) 

Execution 
Time  (Sec) 

FFTIl 

47 

ICO 

2:1 

2.27 

FFTI2 

36 

161 

2N 

2.5 

FFTI3 

36 

112 

2:1 

2.185 

FFTI4 

42^ 

1543 

2N  + N/4  +12 

1.95 

FFTI5 

31^ 

1263 

2N  + N/2  + 12 

) . 20 

^Subroutine  COSINE  has 

10  statements. 

29  words 

^Subroutine  COSINE  has 

8 statements,  19  words 

^Subroutine  UNSCSl  has 

29  statements. 

125  memory  words 

NOTE: 

Subroutine  UNSCI 
(Appendix  F) . 

11  is  faster  than  Subroutine  UIISCR 
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Figure  42 » Subroutine  FFT15;  Threshold  - 90; 
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VI.  Conclusions  and  Recommendations 


Conclusions 

The  major  conclusion  of  this  thesis  Is  that  no  other  orthogonal 
transform  should  be  substituted  for  the  Fourier  transform  in  the 
application  of  digital  signal  processing  techniques  in  TTI  radar 
Imaging.  The  Karhunen-Loeve  transform  has  no  fast  computational 
algorithm.  The  Hankel  transform  was  not  usable  for  Implementation 
since  It  Is  a two-dimensional  transform.  The  Mellln  and  Cosine  (sine) 
transforms  employ  the  FFT  to  compute  the  transform,  thus  no  speed 
advantage  would  be  realized  for  either  of  these  two  transforms.  It 
was  found  that  the  Walsh  power  spectrum  computed  using  the  fast 
Walsh  (Hadamard)  transform,  did  not  Isolate  the  Individual  scatterers 
of  a complex  target  as  the  Fourier  spectrum  was  able  to  do.  The  con- 
version from  the  Walsh  sequency  domain  to  the  Fourier  frequency  domain 
was  found  to  be  computationally  excessive  and  more  than  offset  the  high 
speed  advantag-*  of  the  FWT/FHT.  The  assertion  made  by  Robinson  that 
there  exists  a linear  transformation  between  the  Walsh  power  spectrum 
and  the  Fourier  power  spectrum  was  found  and  proved  to  be  Incorrect. 

The  fixed-point  FFT  algorithm  was  the  fastest  Implementation  for 
TTI  signal  processing.  It  was  considered  that  the  Increase  in  speed 
without  a significant  Increase  in  error  was  significant.  It  is 
believed  that  execution  will  be  even  faster  when  programmed  on  the 
DEC  PDF- 11/40  minicomputer. 

Racommendatlone 

It  Is  recommended  that  the  scale-lnvarlance  property  of  the 
Mellln  transform  using  exponential  sampling  and  the  FFT  be  studied  to 
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determine  If  it  is  feasible  for  compensating  for  the  nonllnearltles 
of  the  Doppler  shift  frequency. 

It  Is  recommended  that  hardware  Implementations  of  the  Cosine 
(Sine)  transform  be  Investigated  for  Implementation  Into  the  TTI 
Imaging  system. 

. Although  no  specific  system  has  yet  been  designed  and  built,  some 
theoretical  systems  have  been  designed.  They  utilize  very  sophisticated 
tactical  radars,  but  well  within  the  state-of-the-art.  It  Is  recom- 
mended that  these  or  possibly  other  systems  be  designed  and  evaluated 
identifying  important  parameters  and  operating  characteristics. 
Additionally,  noise  and  statistical  analysis  of  such  systems  should 
be  conducted. 
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Appendix  A 


Main  Progr^im  TGTIP 

Program  TGTID  reads  in  the  radar  parameters  RLOC,  PRI, 

STARTF,  and  DELTAF  which  describe  the  location  of  the  radar  relative 
to  the  target  and  the  waveform  desired.  NBURST  and  NFPJiQ  are  the  number 
of  bursts  transmitted  to  the  target  and  the  number  of  frequencies  with 
each  burst,  respectively.  Two  other  parameter?  HSCAT  and  SCOR!-!  describe 
the  target,  where  NSCAT  is  the  number  of  scatterers  within  the  target 
and  SCPRM  is  the  set  of  seven  parameters  of  each  scatterer.  The  first 
three  parameters  are  the  coordinates  of  the  scatterer  on  the  coordinate 
system,  o^  Is  the  fourth  parameter,  and  the  last  three  are  cosine 
weighting  terms. 

Program  TGTID  calls  Subroutine  CINIT  which  computes  eight  additional 
parameters  based  on  the  input  parameters.  They  determine  the  size  as  a 
power  of  two  and  other  characteristics  of  the  time  and  frequency  arrays. 
Subroutine  SLANTV  is  called  which  returns  the  simulated  target  Doppler 
returns.  Then,  subroutine  CIMAoE  processes  the  Doppler  signal  returns 
and  displays  the  synthetic  image  by  calling  either  of  the  display 
subroutines  (Appendix  G) . 

Program  TGTID  finally  prints  out  the  radar  and  scatterer  parameters 
along  with  three  calculated  values;  Number  of  Words/Record,  Number  of 
Records,  and  the  maximum  power. 
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»r5  3rt7  = Tfjouj,XAP"r  = 3UTP!J'^,raf»C:«>rT4PE8,  PLOJl 

«JTM'JLAT£  RAIAt  IMAGE  o^T'IPM 

~01Mnn  /^atAR/  PL0':(3)  ,PPI,S‘^A?TF,nELTAF 
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CG’MTJ  /eA=>AM/  oi»C,Hl,H2 
-ommon  /units/  NVLT 
GOM“''N  /ooINTS/  NIMG 
GO-'M-'N  DATA(lOaOO) 

'’A'A  D:/r.i:»15923/tn/2997a253  3.  /,  HI/,  J,  .',/,H2/.46/ 
MATA  MV’_T/6/ 
data  HCMG/8/ 

DA’A  NnS7/13000/ 


INPUT  PARAMETERS 
simulate  radar  RETURN 
MUTPUT  IMAGE  PARAMETERS 

READ',  RLOC,PRl,STARTF,OELTAF 
RPAD',  A0,R0,G0,A  r33T,RDOT,GnOT. 

READ*',  NMURSTfNERET 
READ-,  NSCAT 

PEADfE,*)  ( (SCPRMdfK)  ,1=1,7)  ,K=1,NSCAT) 

CALL  CINIT 
CALL  SLANTV 

CALL  C’MAGE(DATA,ND1,ND2,N0T,NTM) 

WRTT'-('’,12)  RLOC,PRI,STaRTF,  CELTAF, 

* A0,B:  ,GC  ,A'>OT,RDOT,GDOT, 

* N3URST,NFRFO, 

* NSPAT,  ( (SC'*RM(I,K)  ,I=1,7),K=1,NSCAT) 
WRITr(7,i3)  RMAX, NWRnS,NRECS 


12  rOPMAT(///,T3a,''’’MAGE  PaRAMETERS",/,T3(J,15  (1H-)  , //, 

* no, ’‘RADAR  PARAMETERS  ="  , 3E10 . 0 ,F13  .6,2E12 . 6,/, 

» -^10, "ANGLE  DATA  =”,SE10.2,/, 

* 'T10,"NO.  OF  SURSTG  =",I10,/,. 

* OF  FRETJENCIES  =”,110,/, 

* T10,“NO.  OF  SCATTFRFPS  =",110,/, 

* T1  0,"SCATTFREP  PARAMETERS  =" , ( T32 , 7FJ.0 . 3)  ) 

13  PORMAT  (///, TIC, "MAXIMUM  POWE®  =",.F10.2,/, 

* TIC, "NO.  OF  WOROS/RECORn  =",IiO,/, 

* T10,"ND.  OF  RECORDS  =”,110,/) 

STOP  t END 
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Appendix  B 


Support  Scbroutlnes 


Subroutine  CINIT 


Subroutine  CINIT  ■jalculates  NBORD,  NBSZ,  NFSZ,  NDSZ,  NDl,  ND2, 
ND3,  and  NIM  based  on  the  valixes  of  NBURST  and  NFKEQ.  NFSZ  is  the 
size  of  the  frequency  response  data  array  and  NBSZ  is  the  size  of  the 
time  response  data  array.  They  are  related  to  NBORD  by  the  following 
relationships 
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Subroutines  SLANT'/,  CTRAN,  and  DOTP 


Subroutine  SLANTV  Is  the  driving  subroutine  that  generates  the 
simulated  "slant  voltages",  that  Is,  the  Doppler  signal  received  by  the 
radar  set.  Subroutine  SLANTV  calls  Subroutines  HAMWGT,  CTKAN,  DOTF,  FFT6 
(Appendix  C),  and  ROLL.  The  slant  voltage  values  are  written  to  a 
temporary  file  (NVLT)  which  Is  returned  to  Program  TGTID  for  processing. 
Subroutine  SLANTV  uses  the  radar  and  target  parameters  read  in  by 
Program  TGTID  and  computed  by  Subroutine  CINIT.  No  other  input  is 
required. 
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CT 
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-FTt{-‘jrn3n,^OLT{l.l)  ,V0LT(1,2)  ) 

L ^r^LL  vni-i 

:•^F(^;VLT)  ( CJOL' ('♦'<’)♦  T = 1,NPS7)  ,<=1,2) 

3r-,iT»40  K}»/lT 

O^'UR^I 


ooooo  O oooo  nnnnnn 


P,N‘=’C) 

or-L  SMO  ,OMT (7) ^r-(T) ,P(3) 


^noPPruSTP  •'?,^S'J';Pn"‘»^ATI0^4 


no  TO  (3D,43)  ♦N'^C 


iNlTIiLI'^':  RHO 

30  MP''=2 

(PMTd  ) ) 
r^rSIN (PHJ(P) ) 

‘f';=SIN(D^T(3)) 

rA-rn?  (phkd) 

PPrCOS (DMT(?)) 

(PHI(3) ) 

RKT(iii)  = r4*r'^  - <;fl*~p»SG 

pn'^dj?)  = 

= <'G»CA  f >A»~P*CG 
1)  = ‘?a*rG  *■  '^A*'*P*SG 

9Hn(2,2)  = CA^CR 

nM'>(2,3)  s ^A'RG  - OA*GO»CG 

OM')(3,l)  = -SG*CR 

OM'^(3,2)  s S" 

PHT(3,3>  = GR^CG 


'lATRTy  MiJLTTDLIGA-^TON 
30MOIITF  M^w  CODRDT'iATES 

40  no  44  T=lf3 
o(T)=C  . 
no  44  K=l,3 

0(r)  s p(i)  4.  PT(K)  •RMOrT,K) 


pcruRN 

Etn 


o r>  o o 


S'JnoniriKE  OOTP  (V1,V2,VAL) 
’EflL  V1(3),V2(3) 


COMPUTE  DOT  PRODUCT  OP  ? VECTORS 

VAL  = C1. 

A=0. 

•isO. 

00  A5  1=1,3 

A=A+V1(I)»*2 

8=p+V2(I)**2 

40  VAL  = VAL  ♦ V1(I»*V2(II 
VAL=7AL/SQRT{A*9) 

PETUPN 

END 
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Subroutines  ROLL,  WROLL.  AND  IROLL 


Subroutine  ROLL  rearranges  the  FFT  exit  array  Into  normal  viewing 
in  Subroutines  SLANTV  and  CIMAGE.  llils  caii  be  done  either  before  or 
after  the  power  spectrum  is  computed.  Subroutines  WROLL  and  IROLL  perform 
the  same  function  for  the  Fast  Walsh/Hadamard  transforms  (Appendix  D) 
and  the  fixed-point  FFT's  (Appendix  E),  respectively,  in  the  appropriate 
CIMAGE  subroutine. 


! 

I 

! 


116 


o o o o o 


njisojTjfjp  POLL  (*^S7,V) 
PT -L  \/(?‘>P.,2) 

C 

C 

C POLL  thf  T*<aC>p 

C 

NS"'2='<';7/? 
no  U4  KrljNF'” 

L = ''+N^7^* 
m 44  I=i,2 
T=  /(<,  I) 

=\/(L,r) 

44  V(l,I)=T 

c 

pc-rupM 


P'nPOJTINE  MP0LL(M«;7,\/) 

c 


POLL  the  IMftGE 


NO-'SsNS?/? 
no  10  k=i,ne72 
L=<*-N';?2 
T=/(K) 

V(‘')=7{L} 
V(L)=T 
C 

10  CO'JTIMUE 
RP7URM 
ENT 
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O O O O 


*?'P?0'JTIN^  IRT-L  (M*;'’,  rv^,  I\/I) 

T'^'’ ( ?~n  5 , TVT  ( ?-■'  ),T1,T2 


»OLL  "n-  IMA'^F 

no  i<=i,m<577 

l=-'+n';'»2 

*1=177 (<) 
tV-'CK)  =I70(L) 
T7-(L)=T1 
-?=!VT (<) 

7\f’^{<)  =I7I(L1 
TV'(L)  =T2 
4t  rO‘JTISUE 

c 

RF*U7'4 

FN? 


( 


Subroutine  CIMGE 


Subroutine  CIMAGE  constructs  the  Image  from  the  simulated  Doppler 
returns  generated  by  Subroutine  SLANl'V . The  data  are  first  Hamming 
weighted  as  the  sample  values  are  read  from  the  temporary  data  file 
(NVLT).  The  orthogonal  transform  is  taken  along  constant  range,  the 
power  spectrum  computed,  and  the  data  are  "flipped"  by  Subroutine  ROLL. 
The  "image"  data  are  then  written  onto  a second  temporary  file  (NIMG). 
Subroutines  PLOTID  or  BUFOUT  (Appendix  G)  are  then  called  to  display 
the  image.  Subroutine  CIMAGE  is  written  in  three  versions;  the  first 
is  used  for  floating  point  FFT  subroutines,  the  second  is  for  Walsh 
transform  subroutines,  and  the  third  is  for  fixed-point  FFT  subroutines. 
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p 


/ 


« 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


';m >»n.iTTMr  rf'iAC-r  nT'^) 

r-*  /LI'^TT/  VTi|3^T,N-‘n';0,‘J9S2,N5- iET, ‘JP0?T,NC-5T 

ro  I'lO'l  - nv,*^w?f 5,Mi<r05 

rT'^ON  /'Jf.'TTS'  “Jv/L" 

f'O-iMaN  /O'^T'IT';/  MIMC 

®=:''L  (/0LTf?rn,?) 


'OM^TP'.irr  TW;  TMCG^" 


roHP'»Tc  WEIGHT  \/'‘GTO’. 

tkjoiit  volts  G!^*- 

PTVOT  _ p,LflM*  "ANGz  TO  OflNE'’ 

FyTE'O  tmp  ftppav 

TAKE  rcT  aLOM'  CONSTANT  RANGE 

CMAP  HAH/"? 

CT-'^'JTE  poW-?  '“PpCTRUM 
cnHP'iTP  cower 

rOMPij-rr  rcT  '•X"''')TION  Til" 

O'lTPUT  THp 

pall  HAMWG'^{WT1  ,W) 

t*tme=C.O 

no  49  M=1  ,NTm 
NCsNOT* (H-i) 

no  42  N1  = 1,>J01 

or.r5{sjvLT)  ( (VOLTd,  J)  ,’'=1,NFS7)  ,Jsl,2) 
no  41  N7=!,Nn3 
o'''’(Nl  ,1,H3)  =V0LT(NGtNT,1) 
i,l  o"-’(Nl  ,?,NT)  =70LT  {M5+N3,?) 

4?  rOMTINUE 

00  46  N’=l,Nn3 
no  43  IrNPijNOS? 

VOLT(I,1)=0. 

4T  VOLTd  ,2)  =0. 

no  44  1=1, NOl  1 

WOLTd,l)=P:n(T,l,N7)*W{T)  ' 

44  VOLTd  ,2)=nCRd  ,2,N7)*W(I)  j 

F-'tmPzC-COMOCC'’) 

INSERT  FA'='T  hepe 

CALL  -‘■T6{NnO'»n^vnLT(l,l)  ,V0LT(1,2))  i 


F’’TMPr«;rC0N0(CO) 
p-TM-=ETIMr-CTIMF 
TTTME=TTIMr+PTiME 
call  roll (M9S7, VOLT) 
on  41  I=1,NRST 

VOLTd  ,1)='?0P^  (VOL^d,!)  *-2+VOLT(I,’)**2) 
4F  oMaxravtflXKo^iax ,voLT(T,i) ) 

WRTTE(NIMG)  (VOL^d,l),  T=1,N3S7) 

46  CONTINUE 
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o r»  o <r» 


t«  P=:MIf4D  *i\JlT 
P- -JIN')  Ninr, 

C 

W3’’Tr{7,iR) 

IF  ro->«i’' (///,iy,"Fe-S="} 

V''’’‘TF(r,  1C)  TTI*^'" 

10  ro JMAT  {iy,*’T«?a‘4F!^0=M  Fy'TUTION  TIMP  I5",F'^.  + ,**  ?F^. 

— TM'^FRT  IMflGF  ''1‘JnLflY  F'J''®OUTIHF  CAL-  (?)  HF?F 

CALL  OLOTin(MTMf;,f4'^7,MTtf,‘J3??) 

CALL  '»irOMT(‘4THrM‘Jn?,flT‘SN8S7) 


N=*CS=NF?*^ 

MM^nS=N??Z 

PFTMRnI 

F^l^ 


Version  1 


oooooooooooooooo 


S«nROUTINE  CIMftGE  (OCP,M01  ,N02,ND3,  MTM) 

COMMON  /LIMIT/  N3URST,N'<OPD,NBSZ,NFPEa,NFORD,NFSI 

COMMON  /IMAGE/  RMAXtNHPOS, NRECS 

COMMON  /UNITS/  NVLT 

COMMON  /POINTS/  NIM5 

P&AL  P:R(N01tN02,N03) 

REAL  VDLT(256,2) ,H(256) 
complex  V0LTC(256) 

real  VDLTR(256) ,V0LTH(256)  ,T(256) ,R'10) 


CONSTRUCT  THE  IMAGE 


COMPUTE  HAMMING  HEIGHT  VECTOR 
INPUT  VOLTAGES 

PIVOT  - SLANT  RANGE  TO  CROSS  RANGE 
EXTEND  THE  ARRAY 

TAKE  FHT/FHT  ALONG  CONSTANT  RANGE 
SHAP  HALVES 

COMPUTE  POWER  SPECTRUM 
COMPUTE  MAXIMUM  POWER 
COMPUTE  FHT/FHT  EXECUTION  TIME 
OUTPUT  THE  IMAGE 

NS5Z1=NBS2»1 

NBSZZsMBSZ^Z 

NBSZM1=NBSZ-1 

NBOR01sALOG(FLOAT(NBSZ2n/ALOGC2.)*.l 
CALL  HAMHGT(N01,H) 

TTIMEsfl.O 
no  48  M=1,NTM 
NS=N03*(M-1) 

00  42  Nl=lfNDl 

PFAO(NVLT)  ( (VOLTd,  Jl  ,I=ltNFSZ)  t 
00  41  N3=1,N03 
P'*"CNlflfN3)  =V0LT(NS*N3,11 

41  P'JR(H1,2,N3)  =V0LT(NS»N3f  2) 

42  CONTINUE 
00  46  N3=1,N03 
00  43  IsNDl,NBSZ 
VOLT(Itl)=0. 

43  VOLT{I»2)=0. 

DO  44  1=1, NOl 

V0LT(I,1)=PCR(I,1,N3)*H{I) 
h4  V0LT(I,2)=PCR(I,2,N3>*W(I) 

DO  50  I=1,NBSZ 

V0LTC(I)=CMPLX(V0LT(I,1) ,V0LTCI,2)) 
V0LTR(I)=CABS(V0LTC(I)» 

50  CONTINUE 

DO  55  I=NBSZ1,N8SZ2 
VOLTRCI)=0.0 
55  CONTINUE 
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c 


STIHE=SECOND(CP) 

I»^SERT  PAST  transform  CALL  HERE 
call  FHTl (V0LTR5NgS7?,N3ORDl,l) 


ETrME=SECOND(CP) 

PTIME=ETIMfr-STIME 

ttime=ttihe*rtime 

00  49  I=1,NBS2 
IPd.EQ.DGO  TO  51 
IF(I.Ea.N8SZ)G0  TO  47 
V0LTR(I)=V0LTR(2*II**2*V0LTR{2*ifl)**2 
GO  TO  45 

51  V0LTR(1»=V0LTR(1)**2 
GO  TO  45 

47  V0LTR(NBS7)=V0LTR(N9S72)**2 
45  RMAX=AMAXl(RHAX,VOLTRCin 
49  CONTINUE 

CALL  HR0LL(N9SZ,V0LrRI 
HRrTE(MlMG)  (VOLTPd)#  I=1,NBS2) 

^6  CONTINUE  i 

48  PEHINO  NVLT 
REMIND  NIMG 

C ■ 

WPITE(7,15) 

15  FORMAT (///, IX, -FHTl-J 

MRITE(Z,10)  TTIME  ! 

10  EORHATdX, "TRANSFORM  EXECUTION  TIME  IS",F7.4,"  SEC.*) 
C 

c---  — insert  image  display  SUROOUTINE  CALL(S)  here— 

CALL  BUE0UT{NIMG,N03,NTH,NBSZ) 

CALL  PL0TI0<NIMG,N03,NTM,MBSZ) 


NPECS=NESZ 

NMROS=N0SZ 

RETURN 

ENO 


Version  2. 
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oorio-^ooooooooooo 


T 


«f"'’30;TI»j--  ‘.T‘nr,r 

ro  'M0>|  ‘.FOcCIjMC*?*’ 

ro  *M0‘J  /iMjf.r/  ^•4;y,‘|^'-''^,NPc-cs 
rf)->yO-]  /iimtts/  ‘JV/LT 
/OOT'JX'T/  ».'TM^ 

®-*L  '’C?  ( ►'Hi  5 

o":-i,  s/OL'^  {’?«.,?) , tj{?r5) 

TMXri^rp  T VOL  T ’(?“%),  T V?',  TT  ( 25  6) 

I»J‘^Er,ER  Tr(2P?>', 


•'M'lT  PJCT  THE 


rn’^PU'^’E  wE;if;HT 

jHOiir  701 

PT^o’'  - PLA'JT  pfl^rt-  TO  CRo^s 

pyTEMn  THP  AopiY 

Tflv'E  TMTEEF^  -HONG  COMSTA  IT  RANGE 
G^/Ao  HALVE'^ 

rOHPilT"  'IAYTm')'^  D'^gro 

rn^o.iTE  Tsrrr,rD  rpT  EXECUTIO'J  '’I'IP 

IJJTOIIT  THp 

raUL  :0«^TMF(TC,NRG7) 

TALL  -Zfl'IWGT  (KHl  ,W) 

•^"TMFsC.O 
cA''TOR  = 2**15-1 
m 48  HrljNTH 

NGsNOX* (M-i) 

nn  42  Nt=l,Nni 

»E2n(sVLT)  ( (VOITC,  J)  ,’'  = 1,NFSZ)  ,J=1,2) 
no  41  N?  = l,ND'» 

PC'’{N1,1,NT)  =7nLT(N5**»-»,l) 

*♦1  PC='(N1  ,?,M7)  =70L’(N'^«-N3,?) 

42  rONTINIIE  

nn  46  N3=1,N0T 
no  43  I = N01,Mn<?7 

vniTd  ,i)=^0. 

47  Vr)«.T(T,2)=C. 
no  44  1 = 1, Nil 

VO'_T(r  ,1)  =oCR(T,l,N7)  •«(!) 

44  VOLTd  ,2)  =PCP  (I  ,?,N3) ‘Wf  T) 
no  50  I=l,NnE7 

IVOLTP (I)=V0LT (1,1) »FAn’0R/N3S7 
T'OLTI  (I)  =V3LT(T,2)  *«^ACT0R/N3S7 
50  rOMTINUE 

GTTHE=SEC0Nn(0O) 
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o o no 


T3.fl'!crT>w)  CALL  ---------- 

C - 

~T  ’'M"=  '■TIM'^-CTTHr 
•'”■»^<^5-TT^CfP7T‘^c 

ra<:t  Tort.t.  ,T'70f.T?,TVCLTi » 

n't  Un  T = l,Mn«;-» 

ILTR  (I)  =‘;^^ie-T{TV'’iL‘^- (T)  **24-IVCI  TI  (!)••?, -15) 

VTLTd  ,1>  =’'\/nLT:?(I) 

4«;  0'«*y=iM4ri  (^M^y  ,ynt.T(T) » 

(Itl) ♦ I=l,'nsZ) 

46  CT'TI'JU' 
a#  O'-'IISD  NVLT 
C 

W5’'T=-C»,15) 

ir  (///,iy,"pr-i3**) 

M'>’T*(7,10)  TTIMr 

1?  po-MftT (iy,**T°ft‘i'>‘^05'‘^  ry'ruTioN  timr  i~”,p'’.4,** 

- — T.j^r^r  iMAr:-  •'•|-'e“OUTINE  CAL- (5)  M"®'! 

Cl'.L  3L0Tn(fJ:M';,».pT,‘ir>%.j3C2) 

'‘A’.L  ^»ir'3(iT(NlMG,MC3,MT  S‘J0SZ) 


NR'^CSsNFS? 

NM-’OSsNR!?? 

“'•TI.IRS 

F*n 


Version  3. 
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Subroutine  HAMNGT 


This  subroutine  generates  a Hamming  weight  vector  (Equati  on  (23)) 
for  sidelobe  control.  Subroutine  HAMWGT  is  used  in  Subroutines 
SLANTV  and  CmCE. 
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r»  o o o 


r 


ro  iMn>j  /3ar>'\M/ 

Wd) 


■XOMnif'*  tipxRHT  vrcTOO 

=»ELl.  rttryc  r^D  «;in'-LORE  REni'CTIOM 

Nn'>=»4/? 

“=‘*024-l 
^=’.*=’T/(M-1) 

W(«)=l. 

4p  1=1, NT? 
y r HI  4.  W3»rT^(3«J) 

»(''♦!)  =y 
40  W(  >-i)  =x 
C 

P=‘T1JRN 

FNT 
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Appendix  C 


Floating  Point  FFT  Subroutines 

Subroutines  SLANTV  and  CIMAGE  (Version  1)  call  an  FFT  subroutine. 
The  data  must  be  placed  Into  the  proper  form  before  the  specific  FFT 
subroutine  Is  called.  Each  subroutine  Is  documented  Internally  and 
will  not  be  discussed  further. 

Reordering  subroutines  are  listed  In  Appendix  G. 
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/ 


ooonoooo 


c<|oc3jTiMr  f 


A = T*lpi|T  AonAY  ic-  crMr>^c«5 
M - j M sfjUMPr^  or  fJAKPLES 
M s L0n?(H) 

TS  = -1,  T^<J^J~P0R^! 

■ S +11  TNyEO'iE  (IJNN0'’'<4LT ZED) 

TI'^TN'i*  .:3 

nTMp»j5iOKj  8(N) 

ro*PL“y  A,I.I,W,T 
HATA 

Ml/  •»s'^/  2 

J=1 

no  30  1*1, NMl 
IPJI.5E.J)  r,0  TO  1C 
Ts  ^ ( J) 

AfDsfld) 
fl(T)sT 
1C  K='IVZ 

?0  T‘='(<.3E.J)  50  TO  ^0 
J=  l-K 
•fsK/S 
00  TO  20 
10  .Js  l*< 

no  80  t=l,M 

LElsLE/2 
tla<'M3.X(l,  0,  0.) 

I^dS)  40, EC, 50 

40  W=''MOlX(C05(pt/lF1)  ,-5T*j(di/lE1)  ) 

TO  6^ 

50  W=''MPLX(C0S(PI/LFl),5IVfPT/LEl)) 

60  no  80  J=l,LP’l 
no  73  I=J,N,L‘" 

10=1 fL El 

fl<TP)=A(I)-T  

70  A(T)sA  (D  d 
80  U='J*M 

RETURX 


oo  oooooort 


-•4Tt  cpr  <?!|n.''.r)!,i7TN“  T'^  M0')r»='IFD  t/F’SIOM 

n-  F-TIA.  I'  FLlMTSf^'F  Th;;  ^IFST  SET  0*=’ 
rx  mui.TT='LI'^S  TM  "MF  OEFIMaTIDN-T^'- 
■^IMF  1, 

TU7  Mflv/r  rnr  rfiMF  HEAVIN';  AF  I»J 

T.,(r  rc-TiA  ‘^lP?QU''T'iE. 

0'!“''FN' IDS'  A(N) 

m-PLEX  A,II,W,T 

HA’A 

H7E=N/2 

N'1i=N-l 

J=1 

m 30  IsljNMl 
’■Pfl.'F.J)  50  TO  IP 
T=A(J) 

A(  l)=A  (I) 

a(^)=^ 

IP  KsSV? 

3r  TP(K.5F.  J)r,n  TO  7r 
j=  ;-K 
Ks</2 
50  TO  20 
3P  Jvi+K 

no  'K  Ialf^ii2 

Io=l«-l 

T=A(T=>) 

A(TO)=A(I)-T 
/»0  A(T)sA(I)+T 
00  60  L=2tM 

Lp=2*tL 

L"l=LE/3 

(l=''MPLX(1.0,0.0) 

41,«‘2,u2 

i*l  W=0MOLX(0n«;(PT/FL0AT(.Fl) ) ,-STN(PI/-LOAT(LFl))) 
50  TO  U5 

4?  W=''1PLX(005(Pr^f^LOAT  a*^!)  » ,SlN(Pl/='.  OAT  (LED)) 
46  no  SO  J=1,LF1 
m 50  I=J|N,tF 
TP=I4L El 
T=A(IP) *U 
A(TP)=A(I)-T 
50  fl(T)rA{I)*T 
5P  U='I*W 
RETUR'^ 

FNO 
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orjooTooooono 


V 


smoDjTXM?-  "'^•"2  (y,*j':Tai;“,M,siGM) 
n-~TNI  'If  7 fi  = TA?L '•■=■' j 

TMD'JT. 

y(?,ic?-)  » ^fiTa  t*j  colu^'4  i 

0^  ■'•10  WHICH  N 15 

*>1  = 

= -1,  rnpua?'  "»AMSP0R:-1 
= +1.  '^'’•'iNS-ORI 

O’lTsilT; 

y{2*i0’-)  » c-n3>iar''-T9aN5P0RwET  tatc  o'jtpu"^  i-j 

• INVP?«:F-'^^AHSP‘0Rhe1  QATA  TJToiJT  TM 

rr,  -DEry  y (?,  W)  , W 
TNTF3EO  R,  5IG‘^ 

H?=N/? 

'=‘L*N='J 

'■’n  33  J=1,N5TAGF 
»=‘J/ 

M=-N2J 
HT-(2**J)/2 
00  ?3  T=1,MT 
.TM-’J=(I-1)»N?  J 

VOIGN=SIGN 

Jf  '==- LlS?J’‘=HTaM->y5TG)^ 

W=  ■'‘1DLy(cc';(TE’7'’)  iSINCTi^MP)) 
on  20  c.sj.,»j^ 

T'^iq  = ;^+TNj2J 

Tri'iqj- C4.IW2J*2 

TO'njs  TGUOl  + M? J 
Tnin3=i«;ijo+>j9 

yC’,I3MO)=y(i,r«;'ni)  «•  W»X(1,I5U52) 

y(?,i^'jn3)=y (iiiG'ni)  - w*x(i,isub2) 

20  continue 

OO  30  Rsl.N 
30  X(t,9) =X(2,o) 

I'^(5IGW,LT,0  ,)  oftiRM 
no  uo  R=1,N 
40  X(^,R)  =X  (1,R)  />^LTN 

PNO 
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. \ 


no 


ODOOOOOOO 


"^in^OJTiKjr  ‘-rr:,  (y,v,M,N,TS) 


X = ^.TIL  n*  INPUT  HATfl 

Y = rMfir.T’''APV  ''0*lPDMr.jj5  or  oflTA 

N = ?**M  = */“  OATA  ?01fr^ 

I?  = -1*  p 0-3 u 

= *1»  INVF-'''^-  TS-AUCPOPM 
'TJ*^T  TALL  ■*''TT'^(y,'',M,N)  PPL1T^(y,N) 
AND  ^PLDTFfV^N)  TO  ppcRDER 

P'^-^FNPION  V{M),YtN» 

1C  LG=1,M 
I.M''!?*  * (M-LC) 

LTy=2*LHX 

2«:!13P/LTX 
13  1^  = 1, L^y 
Ap:=(LM-1)»5:CL 
F='0S(AR^,) 

<^=-PLDAT(TO)  •'^TNtAP;) 

'’D  10  LT=LIX|N,Liy 
Jl=LI-LTy+LM 
I’=J1*-LNX 

TlrXfJ j) -y( JO) 

T'’=Y(J1)-Y(J2) 

y(  <i)  = y(ji)*y(j2) 

Y(  ll)  = Y(  Jl)  tY  (J?) 
y(  J?)  = C»Tj+S*''’2 
10  V(  i2)  = c»T2-S»Ti 

c 

tail  =»TiTS(y  ,y,m,n) 

c 

RF'IRN 

TMT 


1 

1 


oooonoonorjo 


ns'^n  QM  »?:r:UZ*-JCY  AL'G^T'^HH 

X = CD"'=»L'^X  I'nMT  HATfi  .'■.f'cav 
(=;ES'JLT  T''  ■^'TT'jaMrn  t»,  yj 

fj  r MJ-^nrp  TT  r,J-ra  p(. 

T*^  - {nT'>r'~TnN  O'"  T’SHG'^ORM) 

= “1»  c’GawAsn  T = ANS:‘=‘G=>M 
= ♦•If  TN^‘='-GGf  TDA'I'n^M 


ro  ♦RLEX  X(1C  ?a)  ,!I,W,T 
='^  = 3.1415q?^.5 
GO  ?a  L=i,‘^ 

L''-2**  fM  + i-L) 

L'"1=L'=^/? 

U=(l. 0,0.0 
T'‘fi«:) 

F W=:hplX{OP''(PT'FLOAT(l<-1  > ) ,-Slrj(!^I/"LOaT  (LED)  ) 

00  TO  ■» 

Cy  u=OMOLX(OnO  (PI /FLOAT  (L  ED  ) , SI  M ( PI /^i  0 AT  ( L""!)  ) ) 

7 no  20  J=l,L-l 
no  10  I=J,N,LE 
TP=I+LE1 
T=Y(I)+X(IP) 

X(  TP)  = (X(I) -y  (IP) ) •'J 
10  X(T)=T 
20 

HV2=N/2 

NM1=N-1 

J=1 

no  ?o  1=1, Nil 
TFd.OE.J)  00  TO  25 
TrX(J) 

X{P)=X(T) 

X(T)=T 
25  X=NV2 

2F  T^PK.SE.J)  00  TO  30 
J=  J-K 
K=‘'/2 
00  TO  23 
30  J=  )♦■< 

P.-niPN 

ENO 
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ooooooooooooooo 


nr-Tya'^nM  TM  i-r»r^|irkj~v 

7HT5  CT-J-  <;i|3e>'inTThjr  jz  -*-  MQriFIEO 

O'?  FPTUfl,  IT  'HE  LAST  SET 

PO.Dl_ry  M|ii.Tin|_  jr-  jM  -Jr  DEC  IMATIO'4 -!«- 
P?~Tl-NCY  AL'^ORIT'-im, 

Y r r:)MO|.P/  T‘'=*ll'  TATA  a5>~AY 
(RESULT  TT  ■’ET'J^N'T  IN  X) 

N = NJ'^'^ER  DP  lA'A  ?A'«»'L»‘S 

M = 0'’''E?  OP  <;v<;-'ry 

jr;  - (niap^TTON  OP  To'<N':F0‘?M) 

- -1,  PO^WA"?!  Ti?ANPPneH 
= +-lt  I^YP»SF  TOANSPORM 

On-'OLrx  y (1C2^)  |U,W,T 
P''‘  = 3.1i*59?6=:6T’9 

MM i=y- 1 

OO  ?0 

L-=I*‘ (M+l-L) 

L-l=Lr/? 

U= (1,0  »0. 0) 

T'^(IS>  5«6o 

ty  WsOMPLX(COS(OT/'=‘LnAT(LP1  ))  , -SIN  (PI/- LOAT  (LED  ) ) 
r,0  TO  7 

WsOMOLy(COS(PI'c-LOAr  (LTD  ) ,SIN  (pJ/- . OAT  (LED) ) 

7 no  20  J=1,LE1 
no  10  I=J|N,LE 
ia=I>LEl 
*^sY<i)  ♦y (IP) 

y(TP)  = (X{i)-y  (IP)  )’‘J 
1C  X(T)sT 
20 

no  21  1=1, H, 2 
ia=l*-l 

T=Y(T)  ♦xdP) 

y(TP)=x(i) -X (TP) 

21  X(T)=T 
NV»=N/2 
NMlrN-1 
J=1 

no  3G  1=1, NMl 
TPfI.EP.J)  no  TO  25 
T=y(j) 

X(  f)=X  (I) 

X(T)=T 

25  K=VV2 

26  TP(K,SE.J)  GO  TO  30 
J=  J-K 

K='/2 
no  TO  26 
3C  J=U< 

DCTUfJN 

FNO 
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oooooooooooooo 


'f'nsOJTT.jp 


c 

c 

c 


c 

c 

C 


C 

C 

C 


rn'iLry-T'Ji''''Y  «ITHm  O'"  “Onpipo  jr p-jpm 

^T^JF  r.n*:T*.'F  ~O''»»)TflTI0N 

D£-!-a  a^'^T^DT-'r,  jq  9rv'K>5-  ^JT 

A'»"'C!F?5r«; 

<;T':n  = s'omxF’’  nTor^TinM  toa*j'JP‘0'>M  “LAG 
= -I,  FOR  FO-awSR'^  "FANSFORM 
= ♦!,  -0-’  iNi/ES'?-  "oANSFORV 
T"^MF  = P^L'^A  ’IMF 

^ ~ inr.  0=“  rntisTre*  TPftfj^pQpsj  RLOGK 

= POWP=  0-  ? TO  O^’AIN  NMAX 
NM'X  = L-NG7M  QC  PLOCf  V 

ni  'ENEinKJ  Y(MMAX),0G(?),msk<13) 
rn  <DLPy  y,rxr'3,M0LD,XA 
c'>>ii7AL='NFr  (Oyrp^C'') 

=6 . ? 5 ?1  f F?  C r,*  <; I F,‘j /PL 0 1 ( NMA  X ) 

M-x'd)  =MMfly/2 
no  ic  i=2,NDoiv 
10  '(I)  =M5X(3-1) /2 

MM=nM4X 
‘1M=2 

LO  .'P  OVFR  *iO0M  LAYFRS 

00  50  LaYFR=l,NPOW 

NM=NN/2 

NW=0 

no  4C  1=1, 2 
TT=NN*I 

exes  = G'=-X«(2*oT*rjP' FIGN/NMAX) 

W=«-L0AT(NW)*77 
C5(l)  = CnS(W) 

• C'?(2)  = FrN{M) 

nO'^PUTE  ELFMF'JTS  fo^  nO'H  HALYES  0*^  EACS  RLOOX 
no  20  J=1,NM 

t:=ih-i 

TJ=II-NN 
XA=CX0S*X (TI) 

XfII)=X(IJ)-XA 
20  X(^J)=X(IJ)  ♦XA 
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OOOO  ooooo 


TiMP  UP  CFPTr^;  ny  j 


rnMPJT*" 

'^0  30  L0C-?,*mw 

LL=NM-M«;<  (i_n<') 

T-(LL)  3?,7i*,7: 

70  M-i=LL 
■*?  (LOO 

GO  TO  ‘•0 

ru  M'(=MG<  (LOr+1) 

40  OOMTIm'.IE 
t;C  mm=mh*2 

01  PIMAL  PPfl=<PftGwrMC>4T 
"LSO  hjltT»LY  n^LTi  "^TME 

NWtO 

no  90  T=i,MHAy 
K"^1=NW*1 
Hni_0=y  {‘JWj ) 

50  vcwi)  =y(T)*nTXMr 
62  y(T)=H0L0*nTIMr 

oi|up  IJP  grpTTF  oy  1 
00-‘O|JTF  REV/FPFE  ann?r55 

6!?  no  7C  L0C  = l,MO0W 
LLfNW-MS<<LOC> 

T'^<LL)  72,7^,73 
70  MWrLL 

■^2  N'->=MS'<(L0r)  ♦NM 
GO  TO  90 

74  M'-I=M9<  (LOC  + 1) 

90  r.O'JTI'JIlE  

IP  (SIGN)  ICO, 130, 91 

91  P7F=M.(flx 

95  I=1,NMAX 
95  X('")=y  (D/PTN 
130  PPTU^S 
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ooooooooooooooooo 


PP!?*|W?TOPIWP"W^ 


«:'l'>3r).iT’^N':  (L»X,V) 

ni-?rM~T0'I  V(l),  V(l),  T'^(13),  IU(13) 

‘-'iiivaL- N""  (!•''■.,  T'J(9n»  (ICSfISf3)),  ( rD‘=:fI3('*)  ) . 

^ .f , 

* , (IIS.IGt?)),  (I  JF^TSKl'l))  , 

* (T<Ztr<^(ii) ) t (iLS.TSdZ)) , (:'iF,i';(i3) ) 

rTii^ftL-NC'"  f:4L,I''(l)t,  ( I=^u  , 10 ' 2) ) , {rCL,rj(3n, 

* Ct^^L.T’Ki.))  , (le^LtllKS)  ) , ( If^L,  T'J  ( F.)  i , 

* {Tr,t.,T‘i(“)> , (lHL»ii»(9) ) , (irL*T0(3n, 

* (IJL,IJ(lCn,  (IKLtlUdl))  , (ILL,T'J(12)  ), 

* (I“L,TJ<13) ) 

nn’-fl  ^woai/s  .29 


-PT  ! CAST  c-oijpTt-r  t:^a^jSC0»*1 

1 niM£NSIO»JAL,  ‘‘Txro  ©AOIX 
X = o'rAl  “O'^PONcNTS 
Y = HA-JIMA^Y  C0iO0NFMT9 
2 **  L = '-'O.  OF  ELE1FNT5 

L < : OO-oiftrc;  TNVt^S^ 

MO  FOALIMO  IS  OOME 

FO!»T /AL*’*?';'!  V^PTAOLES  AOF  US^O  ^D*2  S3’*TM5 
AoppoyIMA'^'^  ‘■X'rOTIOV  TIME 

TTM~  s 0 ♦ L ‘ 2»*L  (SE:0NDS) 

^ : 9,t.f-5  FOP  L > C 
0 s 9,7e-5  FOP  L < 0 


LP‘4sIA?S{L) 
Hs?*»L2M 
TPfL»  6,5,9 

6 no  7 Jsl,N 

7 V(  J)=»Y(  J) 

C 

6 T«='  (1.7M-1)  15,15,9 

q Ms  *4 

DO  12  K=2,L2N,2 
m/»=m 

Ms-44/U 

P=TW03I/M4 
DO  12  J=1,M 
tm=^TA=F»(  J-1  ,) 
ri-C03 (THFTA) 
51=STM (THPTA) 

r'’=ci»ci-si*si 

C?r2,*Cl»Sl 

C3=C2*C1"S2*P1 

57=C2*S14-52*C1 

JM  '4sJ-M4 

no  12  I = M4,f»,M4 

Jl=i«^JMM4 

J2=J1*M 

J3=J?fM 

J4=J3fM 


laj 


01  -x(jn  ♦vcj^) 

O'*  ty  ( j 1)  -y  ( j 7) 

r7-Y  ( j ^ ) 4.V  ( j-7) 

o-'^ -Y(Jl)  -Yf  J’) 

o--yf J7) +v( jL) 

o'-.7y  ( j ?)  -y  ( j:,! 

O’  -Y  ( J 2)  4-y(JL) 

05  :Y  ( J’)  -V  ( 

Y f )1 ) = oi  4.cr 
Yf  il)  = ’37fr’ 
oo{J-l)i:,l.l,io 

•■  '4  -'7  = 7 ,•*  - 3 A 

Y ( 13)  = T^oi*  ’ 1 

Vf  17  ) = TH  07»  •■  ;<  _T  Mn^* 

-■■<->1  = 7 l-oq 
-•17?=7y-77 

a(  J7)  = Tm71*~7+T‘;  = 7»<17 
Vf  |7  ) = TM  <7  P'  - 7 -T  '<0  1 * “ 7 
-M  ’'1=7  7_  -p  a 
-‘■<77  = 7u+7f. 

Y(  )4)  = T>i  JI*  ^ 7 4.x,.0  7»  j;i 
V(  iu)  = T'vi3p»';7_T^n^*r;3 

i~n  TO  17 

11  y(  n)  = o24-rs 

Y ( 13)  = 04-0C, 

y ( |7)  = PI -cr 

Y ( |7)  = P3_07 

Y f )4)  =03-05 

Y f )i+)  : D;»fDe, 

13  m ITISIJO 

1'^  T0'L?\’-3*  (L3M'?)  ) 

1P>  <'n  17  t=i,m,3 
oi-y  (I ) 4-y  (T«-l) 
o-5  = y f : ) -y  (T  + i) 

07=Y(I ) ♦ Y (Ifl) 

"’‘♦  = Y(I  )-y  (T  + i) 
v(  ’•)  =^l 

Y{  -)  =03 

Y ( T + l)  =03 

1’  Y { 14-1)  =<7^ 

18  T‘'3=‘)/3 
T”L  = N 
J-13 

m 2':  JJ  = 1,11 
T'3  rj)  = 1 

TOfIS( J*l) -1 .LO.C)  ;o  TT  19 
T'^fj)  = ro(j4-i)  /2 
1«  T'l  (J)  = (j  + 1 ) 

J=  «-l 

30  CO'JTI'JU’I 
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) 


J'^  = o‘ 

n'>  2?  lAri.THL 

?2  I'’  = IA,  T^L  ,‘^'>'5 
22  rr  = I'’ , I^L  ♦ 

on  ??  in  = TC,Tni.,Tn<; 
m P2  :="L,T':~ 

22  = 

•'o  2? 

'>0  22  H = T';,  T‘JL  »T•^•; 

22  TI  = TM,r7j_  ,TTS 

no  9?  TJ  = TT,rJL,T  J<? 
nn  2?  I'^  = I.J,  7<L,TK*; 
nn  22  TLsTk',  Il.L  ♦TL'? 
'^O  22  IH  = TL,  HL  ,IMS 

Trrji-TM)  ->1 

21  t=v(jt) 

V(  IT)  = y<lM) 

X(T‘1)=t 

■^-''(JT) 

V(  IT)  = V(I“) 

Y(TM)=T 

22 

Y'^(L)  7?,7F,3‘^ 
nn  3(b 

jr,  vf  ns-v(  j) 

1"  orrijfjfi 

f^n 
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Appendix  D 


FWT/FHT  Subroutines 

The  fast  Walsh/Hadamard  transforms  listed  In  this  appendix  are 
Internally  documented.  Both  the  input  and  transformed  output  arrays 
are  real  and  in  sequency  order.  The  call  statement  must  list  all  of 
the  parameters  required  by  the  subroutine  parameter  list. 


oooor^oooo-^oooo 


«:m  »5:djti‘IF  "w’l 


'MTl  TRaWSFORH  3° 

IN\/F?5r  f'p  ftS'  IMPUT  tff^TOR 

{\  tra-trt.r''  ^TAL 

n-'SRI  PTIOM  O''  VA'IS^L^'^ 

r - fiPilT  '/FfTno  OP  L®‘»!GTH  N 
N - ‘jijM^rs  OP  para 
M - LO^  PftPP  *’  Qp  *1? 

9*»M  = ‘10.  PP  '^QPPP.  P?R  TH"  I»4TE®VAL 
ry  . POP  -9''UAf»P  TPANPPORHt  X = 1 
- PQR  JNVER5F  T3fl>grP0RMt  X = D 


PT •'EM'; TON  p(Ni 
L=*' 

‘'si 

•'PflX.LT.D  GO  TO  UQ 
m 33  NNsl,^ 

Ts : 

Ls!./2 

no  ?o  NL=1,L 
no  13  M<sl,< 

! = '♦! 

Js-f^K 

P(’)s(r(i)+F(J) )/2. 
10  P(J)sP(I)-P(:) 

2C  T=  I 
30  Ks'/»2 

no  TO  80 

4+0  no  TO  NMsi,M 
T=“ 

L=t./2 

no  60  NL=1*L 
no  50  N^si^K 

T=-r»i 

JsT^X 

F("’)sF  <I)  ♦F(J) 

50  P(  J)=P<I)-P(J)-P(J) 
60  1=1 
70  <=>'♦2 
80  9PTU=’N 
FNn 
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OCJUUOUOUUUCJUUUOOU 


•s’ns^OJTIMf  -M"*  (y,M,M,i<’TPT,MUN) 

T!:^‘'>4"T0MaL  HaT'*M*Pi  tp^NSP'Dv'I 

''DDT  = 1 ! ‘^C  = A'-'^L'^  AMP  Ihl^c-RT 
= ? • TQ  wTTMouT 

'IR  TO  IUV?-~T  WIT-'OUT  3:;‘^aMnLiNG 
= 7 j in"'  UNSCPa E 

NIJM  t NOR-^AL"''^"")  I'NNOPSALITE")  OUTPUT 
= 1 t M^P^ALI■^F'^  (FO^WAPO) 

= 0 * M*JM0RMALI‘"F"»  {INVERSE) 


/API 

X - IM9MT  TA-'S  aoPftX 

M - ‘JIIMTP-R  ;5p  qTTS 
N - NIIM'’F-5  nc*  TATfl  JOINTS  I 


n^'^ENiION  X(‘n 
TPfK0O’’,Uf.l)  GO  TO  71 
HTttsS  , 

irG-7r=4 

no  1 ynlTr^jM 

KH  ILFrKPITC'/? 

KH‘>L‘=‘ls'<HAL''  + l 
00  2 K = '^MftL^l,KST7E,2 
no  2 <1  = <,N,XSI*’E 
nu-tsy/Kl) 

X(<1)=X(<1+1) 

XC^lfDsOUM 
? Cn’ITI'’/UE 
j'f‘^Pr<«:i2F/!t 
v’‘5-’flN=<ni7E 
m 11  NEIT=2,'4PIT 
KG'^l  = J!J*<o*l 

Pn  12  KST=<ST1,  XSIZ'-,X«:pan 

i<GTNMy=<3T+JIIMO-l 

on  12  yEI/'=<FT,  KSINXX 

no  12  XGlrOrK^INfNtXEI'''^ 

n'J-*=X(  KEING) 

XC^GTXG)  rX/KGT'JOJUXO)  • 
y/'^SIMG*  JUMD)  =nUM 
12  Cn‘'TIX!JE 
j:j  •'P=j'i*id/2 
K75AN=<G0AN/2 
11  ro‘'TIMUE 

KGT7E=<ST7P  + K';iTE 
M7:T=MqiT+l 
1 no'JTI*<UE 
71  Cn‘'TI^UF 

mgpan=n 
NnrST=N5PAN/2 
on  ?0  1=1, M 
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?1  L = l,Nni'5T 

21  J=L, 
ni|  -=X(  j) 

X f I)  =TJ'1+V{  j+ijrj 

V(  f+NUPDrnil  l-XfJ+M^T*;*) 
21  rOJTIMI.'F 
M'?''aN=MTIPT 
Nn*ST=N?oaM/? 

2P  rn'iTiviip 

?2  1=1, 

X{^)=y  (D  /(C’LIATCJ) 

2?  '^TJTI'^IUE 

T'^(<0PT,Nr.2)  aE-^IPM 
X‘;T2E=M 
. MT^Ts‘1 
m 3 ><sit=2,h 

<M.1LF=  <SI7p/2 

KWiL^l  =i<HaLf»-l 
J'P'Osl 

no  1?  mi?tt-9,moit 

X-TlsJU'^o  + i 

lu  KSTrXTTl,  KST7=-,XEf'AN 
K'^tN'^X  =<ST+JUM3-l 
'IT  14  KSIMsK'^T,  XPXh4‘<X 
TT  14  KSIMT=t<'?TM,N,<SI7‘* 

nil  i=x{X';iwg) 

y(''SI'4G)rX(<«5lNG  + JUMO) 
X(''SI'4G  + JU‘^o)  aOUH 
li.  CT'lTIMUP 

Ji/  '0=JU'1D+JUMO 
•<’':'’ftN=x«:PANfi<SoA>i 
!■»  ra‘4TINUE 

''0  4 ‘<=KHALri,KGIzr,? 
nn  i*  <l^<t*UK^l7F 
DIJ-*=X(K1) 

X(^l)  = X(<1+1) 

X(<H-1)=0UM 

4 rOMTIMUE  — 

KG''ZE='<MALP 
Mptt='1PIT-1 
3 rOMTiXUE 
PETURN 


OOOOOC-IOOOOOUC^OOO 


^.nOOJ’lNF  r'JTl 


T^vAMFrO^H 


Hl-i  QIJO-53IJT jKir  PX*}T<?  THE  WALS>- 
IN  tme  FOLLOWIN';  WAY 
1.  G'“^JFR5T~'5  A '^FT  D“  WAL^H  FLIN^TIDN*; 

■*.  A3pAr,^jr,F^;  jwru  X'J  a HAOAMARO  'iAToiX 
•'.  'J5FS  TT  TO  rxvjT  T'^F  FORWA^^n 

XM,/E^-r  T3A*'FC’q^i?<  qp  j^p{|t  »/f‘»t'TP 

NO’-E  - msvxmijm  riTTr  Vr.TTP 

■^■FI^jxTIOM  YARTARLFR 
P - IMTUT  v/'^CTOf 
Y - OUTOJT  7CJT03 
N - fj'jMPrp  Qp  narg  F?  moles 
M - OROrp  r)T  '^vjrrwtj  = 

"'X  - FQPWAPO  TP-sj5cnDM,  IX  = 1 
” TMVFOS"  TRAN^oppyi  _ jj 


ni-E^iSIOM  '^{N),Y(N> 
ro  «M0‘J  MAC(1?S123) 
IM’FGEo  WAL 
C 

GAI.L  MALSH(NJ) 

C 

T^flX.FO.C)  00  "^0  ICO 

c 

'•O  ?0  1=1, w 
V(T)=0. 

"G  25  J=1,N 
A=  -»AL(T,  J) 

2S  V( ’■)  s < F(  J)*A  ) +7  (T) 

•»"  V(’’)a\/  (I)  /FLOAT  (N) 

GO  TO  lie 

c 

C iMVEOSt  tPAkjoooom 

IOC  no  73  1=1, M 

^(n=o . 

f'O  70  J=1,N 
fl=''AL(I,  J) 

7C  P(’)=(V(I) »A) ♦r(x) 

Ilf  OPTyi^^j 

ENT 


O O O T 


”^S  5 v--  0“  'JAL3M  FJ^cmsi 

II3VJ3  rn-  M-Tu^n  ->r::--T ',rr,  py  SPJFRaO'J 

rn  's*ovj  WM(1  12F) 

T-!TF'-,-D  t}5L,<3(’) 

LM-INT  ( ( - ( =^L.  n'sT  ( N)  ) / •'  log  (2  . ) ) + .1) 

EG  1=1,^ 

k’='’-i 

J= ; 

f.  Ji=j 
j=  i-t-i 

TT  (1<,  9)  ) 93,03^4,) 

1C  -^-1 
°(  M=-l 
''0  TO  33 
23  Of  J)=l 
3C  K='/2 

T-  ( J.L  E.  1)  '■^n  TO  o 

"f  'D  = ^ ( J)  *0  ( (i  ) 

■"-fL‘4-J)  t.3,..a,0 

4C  WV  (I,  1)  =1 
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Appendix  E 


Fixed-Point  FFT  Subroutines 

The  fixed-point  FFT  subroutines  are  called  by  Subroutine  CIKAGE 
(Version  3).  The  complex  Input  data  are  integer  scaled  values  placed 
into  two  separate  arrays  declared  in  the  calling  program.  All  the 
subroutine  listings  contain  documentation  defining  the  variables 
employed. 

FFTIA  uses  the  values  generated  by  Subroutine  COSINE  to  perform  the 
butterfly  computations.  Subroutine  COSINE  Is  called  once  In  the  calling 
program  and  the  trlgunometric  values  are  passed  as  an  array  In  the  call 
statement  parameter  list. 

Subroutine  FFTI5  uses  the  values  generated  by  Subroutine  COSINE.  It 
is  called  once  In  the  calling  program  and  the  array  passed  through  the 
parameter  list  In  the  call  statement. 

Subroutines  FFTIl,  FFTIA,  and  FFTI5  call  either  Subroutine  UNSCR  or 
Subroutine  UNSCRl  (Appendix  F)  tc  reorder  the  scrambled  Integer  Fourier 
coefficients. 
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Appendix  F 


Reordering  Subroutines 

The  reordering  subroutln  : listed  In  this  appendix  unscramble 
bit-reversed  or  scramble  Into  bit-reversed  order  the  output  sequence  or 
the  Input  sequence,  respectively.  The  exchange  operations  are  performed 
In-place. 

Subroutine  KBITS  Is  written  to  perform  bit  reversal  of  real  arrays. 
It  must  be  called  twlce-once  for  the  real  component  and  once  for  the 
Imaginary  component.  It  can  easily  be  modified  to  handle  complex 
FORTRAN  arrays. 

Subroutines  UNSCR  and  UNSCRl  are  written  to  unscramble  Integer 
Fourier  coefficients.  Subroutine  UNSCRl  Is  modified  from  Subroutine 
KBITS.  Both  subroutines  cap  easily  be  modified  to  unscramble  real, 
complex,  or  two  separate  real  arrays  depending  on  the  need. 
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\ppendix  G 


Display  Subroutines 

Subroutine  PLOTID  creates  a CALCOMP  display  of  the  radar  image 
processed  by  Subroutine  CIMAGE.  It  reads  data  from  the  temporary  image 
file  (NIIIG)  and  compares  each  value  with  a threshold  and  prints  a 
symbol  if  the  valur  exceeds  the  threshold. 

Subroutine  BUFOUT  creates  a line  printer  display  of  the  radar 
Image  generated  by  Subroutine  CIMAGE.  The  procedure  is  the  same  as 
in  Subroutine  PLOTID,  except  the  display  is  printed  on  a hard  copy  by 
an  on-line  printer. 
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Appendix  H 


Radar/Scatterer  Parameters 

The  radar/scatterer  parameters  used  in  this  thesis  are  listed  In 
this  appendix.  Data  Set  1 Is  listed  on  page  and  Data  Set  2 Is 
listed  on  the  following  page. 

They  are  In  the  Format  specified  by  Program  TGTID. 


IMAGE  parameters 


000*0 

000*0 

000*0 

000*0 

o 

UJ 

o o 

«0  o o A0 

e o 

tO  o e iO 

o • 

40  O o CO 

o o 

• • • • 

o 

1 o o 

• 

o 

♦ o 

o o o e 

UJ  O 

o o o o 

o • 

O o lA 

o o 

• • • • 

o 

o 

a* 

a» 

• 

O (Vi 

O O o o 

o o 

o o o o 

o • 

o O o o 

o 

• • • • 

lA 

U\  (VJ  (VI  fO 

o 

• 

* o 

o o o o 

o o 

o o o o 

# 

o o o o 

o 

• • • • 

e o e o 

• o 

O a 


a o o o 
e a o o 
o o o o 

• • • • 

ro  a a ▼< 

I ■ 


• ^ e o ^ o 

0 0^0  o 


in  I 


o o o e 

• • • • 

o in  o a 


II  II  II  II  n H 


t/i 

or 

UJ 

»- 

bJ 

X 


f/l 

o' 

v>  ui 
bJ  int- 
M a u.‘ 
o bj  r 
z or  < 

CT)  UJ  UJ  or 
*-  3 I-  « 
V3  o H-  a 


«r  A'  u.'  ar 

O' »-  3 O'  o or 

« « CD  h.  to  liJ 

CL  O Of 

^ U.  U.  b.  U* 

O'  ui  r>  o o •— 


O U>  • a a « 

Z O O O O 
O'  «t  z r ar  CO 


II  H N 


:>ti.  lb 

zoo 


X a a 

«r  o o 
z r r 


Data  Set  1. 


163 


IMAGE  PARAMEIERS 
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%e  High  Resolution  Radar  Branch  of  the  Rome  Air  Development  Center  has  developet| 
a tactical  tar^t  identification  (TTI)  pulsed-Doppler  radar  system  which  generate: 
two-dimensional  ^Images^  of  aircraft.  The  signal  processing  technique  utilizes 
the  fast  Fourier  transform  (FFT)  to  produce  a slant-range  versus  cross-rang 
display.  If  the  TTI  system  is  to  be  effectively  employed  in  an  aerial  warfare 
environment  then  real-time  processing  is  necessary.  In  an  effort  to  speedup  the 
signal  processing  several  alternative  transforms  v^ere  studied  as  possible 
substitutes  for  the  FFT.  The  Karhunen-Loeve , Cosine  (Sine),  Mellln,  and  Han!:cl 
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^transforms  were  Investigated  and  found  to  be  Infc  uslble  for  use  In  TTI  Imaging. 
The  Walsh  (Hadamard)  transform  was  studied  In  d .all  and  tested  In  a simulation 
program  and  found  that  It  could  not  be  utilized  In  the  TTI  signal  processing. 

IVo  methods  of  converting  from  the  Walsh  sequency  domain  to  the  Fourier 
frequency  domain  were  studied.  The  first  scheme,  a recursive  relationship 
between  the  arithmetic  and  logical  autocorrelation  functions  as  presented  by 
Robinson  was  discovered  to  be  Incorrect.  The  second,  a method  of  computing  the 
Fourier  coefficients  from  the  Walsh  coefficients  of  a function  was  demonstrated 
to  be  too  time  consuming  to  be  Implemented  In  TTI  signal  processing. 

Several  floating-point  FFT  Implementations  were  tested  using  the 
simulation  program.  Also,  several  fixed-point  FFT  algorithms  were  derived  and 
tested.  All  of  these  were  evaluated  on  the  basis  of  speed  and  memory  require- 
ments and  one  fixed-point  FFT  algorithm  was  shown  to  be  fast  enough  and 
accurate  enough  for  Implementation  on  the  TTI  Mln^  puter. 
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